logicaffeine_proof/

lyapunov.rs

//! **Lyapunov-function synthesis** — don't search for the proof, *discover the measure that collapses
//! the problem*, and let the proof fall out as a corollary.
//!
//! This is the thesis of the whole certified-symmetry campaign made into one procedure. A
//! symmetry-breaking refutation is a Lyapunov function for the search dynamics: a scalar potential
//! ("active items remaining") that strictly descends to the goal. Once you have that potential, the
//! refutation, its correctness, and its `O(n²)` complexity certificate are all just *readings* of it
//! ([`crate::complexity`]). So the engine here inverts the usual game:
//!
//! 1. **Synthesize** — search a bounded, polynomial class of candidate potentials. The class is
//!    "covering layouts": a factorization of the variables into `items × bins` under which *swapping
//!    two items is a symmetry of the formula*. Testing a candidate is cheap (one automorphism check),
//!    so the whole search is polynomial — categorically unlike the exponential graph-automorphism
//!    detection it replaces.
//! 2. **Collapse** — if a potential is found, drive the steered Heule descent with the discovered
//!    item-swap symmetry ([`covering_collapse`]); every step is PR-self-checked, fail-closed.
//! 3. **Fall out** — the result is a [`RankedRefutation`] that certifies both correctness *and* its
//!    own polynomial size.
//!
//! When no potential in the class works, the answer is an honest, bounded **impossibility**: "this
//! formula has no covering-symmetry collapse" — `None`, never a wrong answer.

use crate::cdcl::{Lit, SolveResult, Solver};
use crate::complexity::RankedRefutation;
use crate::proof::{Perm, ProofStep, Witness};
use crate::symmetry_detect::{perm_is_automorphism, AutomorphismIndex};
use crate::xorsat::XorEquation;

// =================================================================================================
// The Lyapunov certificate — the unifying object, made rigorous.
// =================================================================================================
//
// A refutation that collapses an exponential search is, formally, a discrete dynamical system whose
// trajectory carries a **Lyapunov function**: a scalar potential `V` over states that is bounded
// below, never increases along the trajectory, strictly decreases across its level set (no infinite
// stall), and reaches its minimum exactly at the goal (⊥). The same four axioms that prove a control
// system converges prove a refutation terminates — and the descent *rate* bounds the proof size, so
// the certificate is simultaneously a termination proof and a complexity bound.
//
// This is the bridge nobody names: termination measures (program verification), ranking functions
// (complexity), and energy/Lyapunov functions (dynamical systems) are the *same* object. Below we
// make it a first-class, machine-checked certificate, and — the load-bearing generalization — we
// show it covers BOTH of our collapse mechanisms: the **geometric** collapse (symmetry: `V` = active
// items remaining) and the **algebraic** collapse (parity: `V` = dimension of the unsolved GF(2)
// system). One framework, two physics.

/// A machine-checked certificate that a refutation's trajectory carries a valid Lyapunov function,
/// with the four dynamical-systems axioms verified explicitly and the resulting complexity bound.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct LyapunovCertificate {
    /// `V` at the start of the trajectory.
    pub initial: u64,
    /// `V` at the goal (its minimum, reached at ⊥).
    pub minimum: u64,
    /// Number of distinct potential values — the length of the strict descent.
    pub levels: u64,
    /// The most work (steps) dissipated at any single potential level.
    pub max_dissipation: u64,
    /// The size bound the descent certifies: `levels · max_dissipation`.
    pub size_bound: u64,
    /// The actual number of trajectory steps (`≤ size_bound`).
    pub total_steps: u64,
    /// Axiom 2: `V` never increases along the trajectory.
    pub monotone: bool,
    /// Axiom 3: `V` strictly decreases across its level set (no level recurs — guaranteed progress).
    pub strict_descent: bool,
    /// Axiom 4: the trajectory reaches the goal (the refutation closes at ⊥).
    pub reaches_goal: bool,
}

/// Verify the **four Lyapunov axioms** over a potential trajectory and return the certificate, or
/// `None` if any axiom fails. (1) bounded below — `u64` is `≥ 0` by type; (2) monotone
/// non-increasing; (3) strict descent across the level set; (4) reaches the minimum at the goal. A
/// trajectory that stalls (a level recurs) or never closes is *not* a Lyapunov function and certifies
/// nothing.
pub fn verify_lyapunov(potential: &[u64], reaches_goal: bool) -> Option<LyapunovCertificate> {
    if potential.is_empty() {
        return None;
    }
    // Axiom 2: monotone non-increasing.
    let monotone = potential.windows(2).all(|w| w[1] <= w[0]);
    if !monotone {
        return None;
    }
    // Axiom 3: the distinct level set (consecutive runs, since non-increasing) strictly decreases —
    // a level, once left, never recurs. This is what forbids an infinite stall.
    let mut levels_seq: Vec<u64> = potential.to_vec();
    levels_seq.dedup();
    let strict_descent = levels_seq.windows(2).all(|w| w[1] < w[0]);
    if !strict_descent {
        return None;
    }
    if !reaches_goal {
        return None;
    }
    // Dissipation per level.
    let mut counts: std::collections::BTreeMap<u64, u64> = std::collections::BTreeMap::new();
    for &p in potential {
        *counts.entry(p).or_insert(0) += 1;
    }
    let levels = counts.len() as u64;
    let max_dissipation = counts.values().copied().max().unwrap_or(0);
    Some(LyapunovCertificate {
        initial: potential[0],
        minimum: *potential.last().unwrap(),
        levels,
        max_dissipation,
        size_bound: levels * max_dissipation,
        total_steps: potential.len() as u64,
        monotone,
        strict_descent,
        reaches_goal,
    })
}

/// Extract and verify the Lyapunov certificate of a **symmetry** refutation: its rank annotation is
/// the potential ("active items remaining"), and the refutation closing is the goal-reaching axiom.
pub fn lyapunov_of_symmetry(ranked: &RankedRefutation) -> Option<LyapunovCertificate> {
    verify_lyapunov(&ranked.ranks, ranked.refuted)
}

/// The Lyapunov function for the **algebraic** collapse — Gaussian elimination over GF(2). The
/// potential is the number of pivots still to be found in the unsolved linear system; each
/// elimination step strictly reduces it, and the trajectory reaches the goal exactly when an
/// inconsistent row `0 = 1` is exposed. Returns the descending potential and whether the goal
/// (contradiction) was reached.
///
/// This is the load-bearing generalization: the *same* Lyapunov machinery that certifies the
/// symmetry collapse certifies the parity collapse — different physics, one potential descending to
/// the goal.
pub fn gaussian_lyapunov(equations: &[XorEquation], num_vars: usize) -> (Vec<u64>, bool) {
    // Each equation is a GF(2) row: a bitset over `num_vars` plus the rhs bit (stored at index nv).
    let words = (num_vars + 1 + 63) / 64;
    let bit = |row: &mut [u64], i: usize| row[i / 64] ^= 1u64 << (i % 64);
    let get = |row: &[u64], i: usize| (row[i / 64] >> (i % 64)) & 1 == 1;
    let mut rows: Vec<Vec<u64>> = equations
        .iter()
        .map(|eq| {
            let mut r = vec![0u64; words];
            for &v in &eq.vars {
                if v < num_vars {
                    bit(&mut r, v);
                }
            }
            if eq.rhs {
                bit(&mut r, num_vars); // the rhs lives at column `num_vars`
            }
            r
        })
        .collect();

    // The potential: number of variable-columns not yet pivoted out, descending as we eliminate.
    let mut trajectory: Vec<u64> = Vec::new();
    let mut remaining = num_vars as u64;
    let mut used = vec![false; rows.len()];

    for col in 0..num_vars {
        // Find an unused row with a 1 in this column — the pivot — and eliminate the column from all
        // other rows.
        let pivot = (0..rows.len()).find(|&r| !used[r] && get(&rows[r], col));
        if let Some(pr) = pivot {
            used[pr] = true;
            let pivot_row = rows[pr].clone();
            for r in 0..rows.len() {
                if r != pr && get(&rows[r], col) {
                    for w in 0..words {
                        rows[r][w] ^= pivot_row[w];
                    }
                }
            }
            remaining -= 1;
            trajectory.push(remaining);
        }
    }
    // After full reduction, an inconsistent system exposes a row that is all-zero in the variables
    // but `1` in the rhs column — the `0 = 1` contradiction. That is the goal; the potential bottoms
    // out there.
    let reached_goal = rows
        .iter()
        .any(|r| (0..num_vars).all(|v| !get(r, v)) && get(r, num_vars));
    if reached_goal {
        trajectory.push(0);
    }
    if trajectory.is_empty() {
        trajectory.push(remaining);
    }
    (trajectory, reached_goal)
}

/// A discovered collapsing potential: the formula's variables read as an `items × bins` grid
/// (`var(i, b) = i*bins + b`) under which item-permutations are symmetries — so the formula is a
/// covering / pigeonhole problem and `items > bins` makes it unsatisfiable.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct CollapsingMeasure {
    pub items: usize,
    pub bins: usize,
}

/// The item-swap permutation over an `items × bins` layout: exchange items `a` and `b` (every bin of
/// one ↔ the same bin of the other). The candidate Lyapunov generator.
fn swap_items(num_vars: usize, bins: usize, a: usize, b: usize) -> Perm {
    Perm::from_images(
        (0..num_vars)
            .map(|idx| {
                let (item, bin) = (idx / bins, idx % bins);
                let ni = if item == a {
                    b
                } else if item == b {
                    a
                } else {
                    item
                };
                Lit::pos((ni * bins + bin) as u32)
            })
            .collect(),
    )
}

/// Drive the steered Heule descent over a discovered `items × bins` covering layout: `bins + 1`
/// items are a tight pigeonhole, so force them out of the bins one at a time, each `¬var(i, bin)`
/// certified by "swap item `i` with the last active item." Returns the ranked refutation; every step
/// is PR-self-checked, so a layout that is not really a covering problem simply fails to refute.
pub fn covering_collapse(
    num_vars: usize,
    formula: &[Vec<Lit>],
    items: usize,
    bins: usize,
) -> RankedRefutation {
    let mut db = formula.to_vec();
    let mut index = AutomorphismIndex::with_clauses(num_vars, formula);
    let mut steps: Vec<ProofStep> = Vec::new();
    let mut ranks: Vec<u64> = Vec::new();
    let var = |i: usize, b: usize| (i * bins + b) as u32;
    let active = bins + 1;

    if active <= items {
        for m in (2..=active).rev() {
            let bin = m - 2;
            let last = m - 1;
            for i in 0..last {
                let clause = vec![Lit::neg(var(i, bin))];
                let witness = Witness::Substitution(swap_items(num_vars, bins, i, last));
                if crate::pr::is_pr_indexed(num_vars, &db, &mut index, &clause, &witness) {
                    db.push(clause.clone());
                    index.insert(clause.clone());
                    steps.push(ProofStep::Pr { clause, witness });
                    ranks.push(m as u64);
                }
            }
        }
    }

    let mut solver = Solver::new(num_vars);
    for c in &db {
        solver.add_clause(c.clone());
    }
    let refuted = match solver.solve() {
        SolveResult::Sat(_) => false,
        SolveResult::Unsat => {
            for lc in solver.learned() {
                steps.push(ProofStep::Rup(lc.lits.clone()));
                ranks.push(1);
            }
            crate::pr::check_pr_refutation_fast(num_vars, formula, &steps)
        }
    };

    RankedRefutation {
        refuted,
        steps,
        ranks,
    }
}

/// Search the bounded class of covering layouts for one whose item-swap is a genuine symmetry — a
/// cheap pre-filter (one automorphism check per candidate factorization). Returns the first such
/// measure in the `items > bins` (unsatisfiable) direction, or `None`.
pub fn synthesize_measure(num_vars: usize, formula: &[Vec<Lit>]) -> Option<CollapsingMeasure> {
    for bins in 1..num_vars {
        if num_vars % bins != 0 {
            continue;
        }
        let items = num_vars / bins;
        if items <= bins {
            continue; // need more items than bins for a pigeonhole contradiction
        }
        if perm_is_automorphism(formula, &swap_items(num_vars, bins, 0, 1)) {
            return Some(CollapsingMeasure { items, bins });
        }
    }
    None
}

/// The whole inversion: **discover the Lyapunov function, then let the proof fall out.** Search the
/// covering-layout class; for each candidate whose item-swap is a symmetry, run the steered collapse;
/// return the first that yields a checking refutation, paired with the measure that drove it. `None`
/// is a bounded impossibility — no covering-symmetry collapse exists for this formula.
pub fn solve_by_measure_synthesis(
    num_vars: usize,
    formula: &[Vec<Lit>],
) -> Option<(CollapsingMeasure, RankedRefutation)> {
    for bins in 1..num_vars {
        if num_vars % bins != 0 {
            continue;
        }
        let items = num_vars / bins;
        if items <= bins {
            continue;
        }
        if !perm_is_automorphism(formula, &swap_items(num_vars, bins, 0, 1)) {
            continue;
        }
        let ranked = covering_collapse(num_vars, formula, items, bins);
        if ranked.refuted {
            return Some((CollapsingMeasure { items, bins }, ranked));
        }
    }
    None
}

// =================================================================================================
// The ⟸ theorem: a poly-description Lyapunov measure ⟹ a poly-size, checkable proof.
// =================================================================================================
//
// THEOREM (⟸). Let `F` be a formula and `M` a Lyapunov measure for it: an initial potential `L`, a
// per-level width `w`, and, from any database `D ⊇ F` at potential level `ℓ > 0`, a set of at most
// `w` clause additions, EACH redundant (RUP/PR/SR) against `D`, that descend to level `ℓ-1`; with
// level `0` forcing ⊥ derivable by unit propagation. Then `F` has a checkable refutation of size
// `≤ L·w + |closure|`, constructible in time `O(L·w·c)` (`c` = per-step check cost).
//
// PROOF (constructive, = `proof_from_measure`). Induct on the potential. Apply the descent steps
// level by level; each adds ≤ `w` certified-redundant clauses and strictly decreases the potential,
// so after ≤ `L` levels the potential is `0` and ⊥ is RUP-derivable. Every added clause is a
// checkable redundancy step, so the concatenation is a checkable refutation, with ≤ `L·w` descent
// steps. ∎
//
// This is what makes the Lyapunov certificate non-trivial: the measure is not *read off* a proof,
// it *generates* one. The hypothesis is the trait `LyapunovMeasure`; the conclusion is the output of
// `proof_from_measure`; the tests instantiate the theorem on structurally different measures and
// machine-check the size bound and the (independent) re-checking of the produced proof.
//
// REVIEWER'S KILLER QUESTION — "is this just resolution width renamed?" NO, and the construction is
// the evidence: the descent steps are PR/SR (substitution-redundant), not resolution steps. So the
// theorem produces a `Θ(n²)` proof of pigeonhole — a formula whose every *resolution* proof is
// `2^Ω(n)` (Haken 1985), at any width. A resolution-width object cannot do that. The measure lives in
// a strictly stronger system; whether its *impossibility* direction beats known lower-bound
// techniques is the open ⟹ converse, stated honestly as such.

use crate::cdcl::Lit as Lit_;

/// The hypothesis of the ⟸ theorem: a Lyapunov measure that *generates* certified descent steps for
/// a formula, independent of any particular family.
pub trait LyapunovMeasure {
    /// The number of variables of the formula.
    fn num_vars(&self) -> usize;
    /// The formula `F` the measure refutes.
    fn formula(&self) -> &[Vec<Lit_>];
    /// The initial potential `L` (the number of descent levels).
    fn initial_potential(&self) -> u64;
    /// The per-level width bound `w` (max certified additions per level).
    fn width(&self) -> u64;
    /// The certified-redundant clause additions descending from `level` to `level-1`, given the
    /// current database. Each MUST be RUP/PR/SR against `db` (the constructor re-checks, fail-closed),
    /// and there must be at most `width()` of them.
    fn descent_step(&self, level: u64, db: &[Vec<Lit_>]) -> Vec<(Vec<Lit_>, Witness)>;
}

/// **The constructive proof of the ⟸ theorem.** Drive the measure's descent from `L` down to `0`,
/// self-checking every step (fail-closed), then close with unit propagation. Returns a ranked
/// refutation whose descent has `≤ L·w` steps — read its certificate with
/// [`crate::complexity::RankedRefutation::certify`].
pub fn proof_from_measure<M: LyapunovMeasure>(measure: &M) -> RankedRefutation {
    let nv = measure.num_vars();
    let formula = measure.formula();
    let mut db: Vec<Vec<Lit_>> = formula.to_vec();
    let mut index = AutomorphismIndex::with_clauses(nv, formula);
    let mut steps: Vec<ProofStep> = Vec::new();
    let mut ranks: Vec<u64> = Vec::new();

    let l = measure.initial_potential();
    for level in (1..=l).rev() {
        for (clause, witness) in measure.descent_step(level, &db) {
            if crate::pr::is_pr_indexed(nv, &db, &mut index, &clause, &witness) {
                db.push(clause.clone());
                index.insert(clause.clone());
                steps.push(ProofStep::Pr { clause, witness });
                ranks.push(level);
            }
        }
    }

    let mut solver = Solver::new(nv);
    for c in &db {
        solver.add_clause(c.clone());
    }
    let refuted = match solver.solve() {
        SolveResult::Sat(_) => false,
        SolveResult::Unsat => {
            for lc in solver.learned() {
                steps.push(ProofStep::Rup(lc.lits.clone()));
                ranks.push(0);
            }
            crate::pr::check_pr_refutation_fast(nv, formula, &steps)
        }
    };
    RankedRefutation {
        refuted,
        steps,
        ranks,
    }
}

/// The **⟹ direction** of the characterization: any checkable refutation of `n_steps` steps induces
/// a Lyapunov measure — rank the steps in descending order. The induced trajectory is a valid
/// Lyapunov function of size `n_steps`, so the minimum measure cost `μ*(F) ≤ (min proof size)`.
///
/// Combined with the ⟸ theorem (`proof_from_measure`: proof size `≤ L·w`), this gives the
/// **CHARACTERIZATION `μ*(F) = Θ(min proof size)`** — the measure framework is *equivalent* to the
/// proof framework, so a lower bound on the measure cost IS a proof-size lower bound. That equivalence
/// is the rigorous foundation of the "no bounded measure ⟹ no short proof" direction (the prize),
/// stated honestly: it makes measure lower bounds and proof lower bounds *the same problem*, rather
/// than giving a new technique to prove either.
pub fn proof_induced_measure(n_steps: usize) -> Vec<u64> {
    if n_steps == 0 {
        return vec![0];
    }
    (0..n_steps).map(|i| (n_steps - i) as u64).collect() // S, S-1, …, 1
}

/// A covering Lyapunov measure: `items` into `bins` (`var(i,b) = i*bins + b`) with item-swap
/// witnesses. ONE measure type that instantiates the theorem for *every* covering family — PHP
/// (`items = n, bins = n-1`), clique-coloring (`items = n, bins = k`), and any other generalized
/// pigeonhole. The "lift": the family is data, the measure is uniform.
#[derive(Clone)]
pub struct CoveringMeasure {
    pub num_vars: usize,
    pub formula: Vec<Vec<Lit_>>,
    pub items: usize,
    pub bins: usize,
}

impl LyapunovMeasure for CoveringMeasure {
    fn num_vars(&self) -> usize {
        self.num_vars
    }
    fn formula(&self) -> &[Vec<Lit_>] {
        &self.formula
    }
    fn initial_potential(&self) -> u64 {
        (self.bins + 1).min(self.items) as u64 // L = active items (a tight pigeonhole over the bins)
    }
    fn width(&self) -> u64 {
        self.items as u64 // ≤ items certified additions per level
    }
    fn descent_step(&self, level: u64, _db: &[Vec<Lit_>]) -> Vec<(Vec<Lit_>, Witness)> {
        let m = level as usize;
        if m < 2 {
            return Vec::new();
        }
        let bin = m - 2;
        let last = m - 1;
        (0..last)
            .map(|i| {
                let clause = vec![Lit_::neg((i * self.bins + bin) as u32)];
                let witness = Witness::Substitution(swap_items(self.num_vars, self.bins, i, last));
                (clause, witness)
            })
            .collect()
    }
}

/// The Lyapunov function for the **cardinality / cutting-planes** collapse — the THIRD physics. The
/// Cook–Coullard–Turán refutation of PHP sums `2n-1` pseudo-Boolean constraints; each addition cancels
/// a literal against its negation and tightens the accumulated constraint toward the infeasible
/// `0 ≥ 1`. The potential is the number of constraints still to combine, descending to `0` exactly as
/// the contradiction is exposed.
///
/// Its significance is **non-uniqueness**: PHP has a covering-symmetry measure (`n` levels × `n`
/// width) *and* this cutting-planes measure (`2n-1` linear steps) — two valid Lyapunov functions for
/// the *same* formula, in two different proof systems. The measure is a property of the
/// *problem-plus-structure*, not a unique object — which is exactly why the framework spans systems.
pub fn cutting_planes_lyapunov(n: usize) -> (Vec<u64>, bool) {
    use crate::pseudo_boolean::PbConstraint;
    if n < 2 {
        return (vec![0], true);
    }
    let holes = n - 1;
    let var = |i: usize, h: usize| i * holes + h;
    let mut constraints: Vec<PbConstraint> = Vec::new();
    for i in 0..n {
        let lits: Vec<(usize, bool)> = (0..holes).map(|h| (var(i, h), true)).collect();
        constraints.push(PbConstraint::clause(&lits)); // Σ_h x_{i,h} ≥ 1
    }
    for h in 0..holes {
        let lits: Vec<(usize, bool)> = (0..n).map(|i| (var(i, h), true)).collect();
        constraints.push(PbConstraint::at_most(&lits, 1)); // Σ_i x_{i,h} ≤ 1
    }
    let total = constraints.len() as u64;
    let mut acc: Option<PbConstraint> = None;
    let mut trajectory: Vec<u64> = Vec::new();
    for (k, c) in constraints.into_iter().enumerate() {
        acc = Some(match acc.take() {
            None => c,
            Some(a) => a.add(&c),
        });
        trajectory.push(total - 1 - k as u64); // potential = constraints remaining to combine
    }
    let reached_goal = acc.map_or(false, |a| a.is_contradiction());
    (trajectory, reached_goal)
}

/// A covering measure restricted to the absolute rounds `[lo, hi]` — it breaks only those covering
/// levels, leaving a residual for a later stage. The tool that makes COMPOSITION non-trivial: a stage
/// that genuinely hands a smaller problem to the next.
#[derive(Clone)]
pub struct PartialCoveringMeasure {
    pub base: CoveringMeasure,
    pub lo: usize,
    pub hi: usize,
}

impl LyapunovMeasure for PartialCoveringMeasure {
    fn num_vars(&self) -> usize {
        self.base.num_vars
    }
    fn formula(&self) -> &[Vec<Lit_>] {
        &self.base.formula
    }
    fn initial_potential(&self) -> u64 {
        (self.hi.saturating_sub(self.lo) + 1) as u64
    }
    fn width(&self) -> u64 {
        self.base.width()
    }
    fn descent_step(&self, local_level: u64, db: &[Vec<Lit_>]) -> Vec<(Vec<Lit_>, Witness)> {
        // local_level 1..=(hi-lo+1) ⇒ absolute round lo + (local_level - 1); as the composer iterates
        // local_level high→low, the absolute rounds descend hi→lo.
        let abs = self.lo + (local_level as usize).saturating_sub(1);
        if abs < self.lo || abs > self.hi {
            return Vec::new();
        }
        self.base.descent_step(abs as u64, db)
    }
}

/// **Compose collapses — the `Poly` wiring diagram on certified descents.** Each stage's descent runs
/// on the residual database of the previous; the potentials are *banded* (earlier stages occupy the
/// higher rank band) so the composite strictly descends across every stage boundary; the whole closes
/// by resolution. The result is ONE refutation carrying ONE combined Lyapunov certificate — two (or
/// more) collapses wired into a single descent. Fail-closed: each step is PR-self-checked.
///
/// Categorically this is sequential composition of coalgebras: stages wire end-to-end, and the
/// banded potential is the single countdown morphism of the composite system.
pub fn compose_collapses(
    num_vars: usize,
    formula: &[Vec<Lit_>],
    stages: &[&dyn LyapunovMeasure],
) -> RankedRefutation {
    let mut db: Vec<Vec<Lit_>> = formula.to_vec();
    let mut index = AutomorphismIndex::with_clauses(num_vars, formula);
    let mut steps: Vec<ProofStep> = Vec::new();
    let mut ranks: Vec<u64> = Vec::new();
    let mut above: u64 = stages.iter().map(|s| s.initial_potential()).sum();

    for stage in stages {
        let l = stage.initial_potential();
        above -= l; // this stage occupies the rank band [above+1 ..= above+l]
        for level in (1..=l).rev() {
            for (clause, witness) in stage.descent_step(level, &db) {
                if crate::pr::is_pr_indexed(num_vars, &db, &mut index, &clause, &witness) {
                    db.push(clause.clone());
                    index.insert(clause.clone());
                    steps.push(ProofStep::Pr { clause, witness });
                    ranks.push(above + level);
                }
            }
        }
    }

    let mut solver = Solver::new(num_vars);
    for c in &db {
        solver.add_clause(c.clone());
    }
    let refuted = match solver.solve() {
        SolveResult::Sat(_) => false,
        SolveResult::Unsat => {
            for lc in solver.learned() {
                steps.push(ProofStep::Rup(lc.lits.clone()));
                ranks.push(0);
            }
            crate::pr::check_pr_refutation_fast(num_vars, formula, &steps)
        }
    };
    RankedRefutation {
        refuted,
        steps,
        ranks,
    }
}

// =================================================================================================
// Stepping past: a UNIFIED auto-collapse engine that recognizes WHICH physics applies.
// =================================================================================================

/// Recover the XOR (parity) constraints latent in a CNF: a constraint `x₁⊕…⊕x_k = b` is encoded by
/// exactly the `2^(k-1)` clauses over `{x₁,…,x_k}` whose negated-literal count has one fixed parity.
/// We group clauses by their variable set and emit an [`XorEquation`] for each group that is exactly
/// such a gadget. Sound: every emitted equation is logically equivalent to a clause group present in
/// the formula, so a refutation of the extracted (sub)system implies the formula is UNSAT.
pub fn extract_xor(num_vars: usize, clauses: &[Vec<Lit_>]) -> Vec<XorEquation> {
    use std::collections::HashMap;
    let mut groups: HashMap<Vec<usize>, Vec<u32>> = HashMap::new(); // var-set → neg-parities
    let mut members: HashMap<Vec<usize>, Vec<Vec<u32>>> = HashMap::new(); // var-set → sign patterns
    for c in clauses {
        let mut vs: Vec<usize> = c.iter().map(|l| l.var() as usize).collect();
        vs.sort_unstable();
        vs.dedup();
        if vs.len() != c.len() || vs.iter().any(|&v| v >= num_vars) {
            continue; // skip tautological / malformed clauses
        }
        let neg = c.iter().filter(|l| !l.is_positive()).count() as u32;
        groups.entry(vs.clone()).or_default().push(neg % 2);
        // a canonical sign pattern (sorted by var) to dedup identical clauses
        let mut sig: Vec<u32> = c
            .iter()
            .map(|l| l.var() * 2 + u32::from(!l.is_positive()))
            .collect();
        sig.sort_unstable();
        members.entry(vs).or_default().push(sig);
    }
    let mut eqs = Vec::new();
    for (vars, parities) in groups {
        let k = vars.len();
        if k == 0 || k > 12 {
            continue;
        }
        let expected = 1usize << (k - 1);
        // exactly 2^(k-1) clauses, all with the same negated-parity, and all distinct.
        if parities.len() != expected || !parities.iter().all(|&p| p == parities[0]) {
            continue;
        }
        let mut sigs = members[&vars].clone();
        sigs.sort();
        sigs.dedup();
        if sigs.len() != expected {
            continue; // not all-distinct ⇒ not a clean XOR gadget
        }
        let b = 1 - parities[0]; // negated-parity p ⇒ XOR rhs = 1 - p
        eqs.push(XorEquation::new(vars, b == 1));
    }
    eqs
}

/// Recover a bipartite covering structure from opaque clauses — the cardinality analogue of
/// [`extract_xor`]. Returns `(rows, columns)` as variable-index groups when the formula decomposes
/// cleanly into "each item is in ≥1 bin" rows (positive disjunctions) and "each bin holds ≤1 item"
/// columns (FULL pairwise-exclusion cliques over a variable set), with every variable in exactly one
/// row and one column. Conservative — a clause that is neither, a variable shared by two rows, an
/// exclusion over an unknown variable, or a column that is not a full clique ⇒ `None`. Unlike
/// [`synthesize_measure`], this needs no `items × bins` factorization and no swap symmetry, so it
/// recognizes *non-uniform / asymmetric* coverings (e.g. the mutilated chessboard).
fn discover_covering(
    num_vars: usize,
    formula: &[Vec<Lit_>],
) -> Option<(Vec<Vec<usize>>, Vec<Vec<usize>>)> {
    let mut rows: Vec<Vec<usize>> = Vec::new();
    let mut excl: Vec<(usize, usize)> = Vec::new();
    for c in formula {
        if c.is_empty() {
            return None;
        }
        if c.iter().all(|l| l.is_positive()) {
            rows.push(c.iter().map(|l| l.var() as usize).collect()); // at-least-one row
        } else if c.len() == 2 && c.iter().all(|l| !l.is_positive()) {
            excl.push((c[0].var() as usize, c[1].var() as usize)); // at-most-one pair
        } else {
            return None;
        }
    }
    if rows.is_empty() {
        return None;
    }
    // Every variable in EXACTLY one row (a clean item partition).
    let mut row_of = vec![usize::MAX; num_vars];
    for (i, r) in rows.iter().enumerate() {
        for &v in r {
            if v >= num_vars || row_of[v] != usize::MAX {
                return None;
            }
            row_of[v] = i;
        }
    }
    // Union-find over exclusion pairs → bin (column) components.
    let mut parent: Vec<usize> = (0..num_vars).collect();
    let find = |parent: &mut Vec<usize>, mut x: usize| {
        while parent[x] != x {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        x
    };
    let mut excl_set: std::collections::HashSet<(usize, usize)> = std::collections::HashSet::new();
    for &(a, b) in &excl {
        if a >= num_vars || b >= num_vars || row_of[a] == usize::MAX || row_of[b] == usize::MAX {
            return None; // exclusion over a variable not in any row
        }
        excl_set.insert((a.min(b), a.max(b)));
        let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
        if ra != rb {
            parent[ra] = rb;
        }
    }
    // Group the row-variables into columns by their union-find root.
    let mut col_id: std::collections::HashMap<usize, usize> = std::collections::HashMap::new();
    let mut columns: Vec<Vec<usize>> = Vec::new();
    for v in 0..num_vars {
        if row_of[v] == usize::MAX {
            continue;
        }
        let root = find(&mut parent, v);
        let id = *col_id.entry(root).or_insert_with(|| {
            columns.push(Vec::new());
            columns.len() - 1
        });
        columns[id].push(v);
    }
    // Each multi-member column must be a FULL clique of exclusions, so "at most one" is genuinely
    // enforced (otherwise the ≤1 PB constraint is not implied by the present clauses).
    for col in &columns {
        for i in 0..col.len() {
            for j in (i + 1)..col.len() {
                if !excl_set.contains(&(col[i].min(col[j]), col[i].max(col[j]))) {
                    return None;
                }
            }
        }
    }
    Some((rows, columns))
}

/// **The third physics, made automatic: cardinality / cutting-planes collapse of a discovered
/// covering.** When the formula is a bipartite covering with more items than bins, summing the
/// "item ≥ 1" rows against the "bin ≤ 1" columns telescopes to `0 ≥ (items − bins) > 0` — the
/// Cook–Coullard–Turán refutation, in the cutting-planes proof system. Sound: every summed
/// constraint is implied by clauses present in the formula (full-clique-checked), so a contradiction
/// proves UNSAT. Returns `(trajectory, reached_goal, constraints)` — the GF(2)-analogue Lyapunov
/// descent (constraints remaining to combine). Crucially this needs **no symmetry**, so it collapses
/// asymmetric coverings the geometric route cannot.
/// Recover the cardinality structure of an opaque CNF as PB constraints — a `≥ 1` clause per at-least-one
/// row and a `≤ 1` per at-most-one clique (the cardinality the pairwise encoding loses). The clause-level
/// recognizer (the [`extract_xor`] analogue for counting) that feeds BOTH the static [`cardinality_collapse`]
/// cut and the live [`crate::pseudo_boolean::CardinalityTheory`]. `None` when the formula is not a clean
/// covering, so callers fall through. Each returned constraint is *implied by* clauses present in the formula
/// (the columns are full-clique-checked by `discover_covering`), so feeding them to a solver is sound.
pub fn recover_cardinality_constraints(
    num_vars: usize,
    formula: &[Vec<Lit_>],
) -> Option<Vec<crate::pseudo_boolean::PbConstraint>> {
    use crate::pseudo_boolean::PbConstraint;
    let (rows, columns) = discover_covering(num_vars, formula)?;
    let mut constraints: Vec<PbConstraint> = Vec::new();
    for r in &rows {
        constraints.push(PbConstraint::clause(
            &r.iter().map(|&v| (v, true)).collect::<Vec<_>>(),
        )); // Σ ≥ 1
    }
    for col in &columns {
        constraints.push(PbConstraint::at_most(
            &col.iter().map(|&v| (v, true)).collect::<Vec<_>>(),
            1,
        )); // Σ ≤ 1
    }
    Some(constraints)
}

pub fn cardinality_collapse(
    num_vars: usize,
    formula: &[Vec<Lit_>],
) -> Option<(Vec<u64>, bool, usize)> {
    use crate::pseudo_boolean::PbConstraint;
    let constraints = recover_cardinality_constraints(num_vars, formula)?;
    let total = constraints.len() as u64;
    let mut acc: Option<PbConstraint> = None;
    let mut trajectory: Vec<u64> = Vec::new();
    for (k, c) in constraints.into_iter().enumerate() {
        acc = Some(match acc.take() {
            None => c,
            Some(a) => a.add(&c),
        });
        trajectory.push(total - 1 - k as u64);
    }
    let reached_goal = acc.map_or(false, |a| a.is_contradiction());
    Some((trajectory, reached_goal, total as usize))
}

/// Recover the at-most-one cardinality **substructure** from arbitrary clauses — the substructure twin of
/// `discover_covering`, which requires the WHOLE formula to be a covering (so a single XOR gadget makes it
/// decline). Binary all-negative clauses `¬a ∨ ¬b` are at-most-one edges; this returns an `at_most(group, 1)`
/// PB constraint per maximal exclusion clique of size ≥ 2, IGNORING every other clause. Each emitted group is
/// a verified FULL clique (every pair has its exclusion clause present), so its `≤ 1` is implied by clauses
/// in the formula — sound to feed a solver alongside a recovered parity system (the par32 fusion).
pub fn recover_at_most_one(
    num_vars: usize,
    clauses: &[Vec<Lit_>],
) -> Vec<crate::pseudo_boolean::PbConstraint> {
    use crate::pseudo_boolean::PbConstraint;
    use std::collections::HashSet;
    let mut adj: Vec<HashSet<usize>> = vec![HashSet::new(); num_vars];
    for c in clauses {
        if c.len() == 2 && c.iter().all(|l| !l.is_positive()) {
            let (a, b) = (c[0].var() as usize, c[1].var() as usize);
            if a != b && a < num_vars && b < num_vars {
                adj[a].insert(b);
                adj[b].insert(a);
            }
        }
    }
    // Greedy maximal-clique cover: grow a clique from each vertex (densest-first), keeping a group only when
    // it covers an exclusion edge no earlier group did. Each grown clique is complete by construction (a
    // vertex joins only if adjacent to ALL current members), so `at_most(group, 1)` is genuinely implied.
    let mut order: Vec<usize> = (0..num_vars).filter(|&v| !adj[v].is_empty()).collect();
    order.sort_by_key(|&v| std::cmp::Reverse(adj[v].len()));
    let mut covered: HashSet<(usize, usize)> = HashSet::new();
    let mut out: Vec<PbConstraint> = Vec::new();
    for &start in &order {
        let mut clique = vec![start];
        let mut cand: Vec<usize> = adj[start].iter().copied().collect();
        cand.sort_by_key(|&v| std::cmp::Reverse(adj[v].len()));
        for &v in &cand {
            if clique.iter().all(|&u| adj[v].contains(&u)) {
                clique.push(v);
            }
        }
        if clique.len() < 2 {
            continue;
        }
        clique.sort_unstable();
        let mut fresh = false;
        for i in 0..clique.len() {
            for j in (i + 1)..clique.len() {
                if covered.insert((clique[i], clique[j])) {
                    fresh = true;
                }
            }
        }
        if fresh {
            out.push(PbConstraint::at_most(
                &clique.iter().map(|&v| (v, true)).collect::<Vec<_>>(),
                1,
            ));
        }
    }
    out
}

/// Enumerate every `width`-subset of `items` (in index order), calling `f`; stop early and return `false`
/// the moment `f` returns `false`, else `true` once all subsets are visited. The shared combinatorial core
/// of [`recover_at_most_k`] — used both to TEST a candidate ("do all `k`-subsets satisfy …?") and to mark
/// the covered subsets of an emitted group.
fn for_each_combo<F: FnMut(&[usize]) -> bool>(
    items: &[usize],
    width: usize,
    start: usize,
    cur: &mut Vec<usize>,
    f: &mut F,
) -> bool {
    if cur.len() == width {
        return f(cur);
    }
    for i in start..items.len() {
        cur.push(items[i]);
        let cont = for_each_combo(items, width, i + 1, cur, f);
        cur.pop();
        if !cont {
            return false;
        }
    }
    true
}

/// Recover the at-most-`k` cardinality substructure over arbitrary LITERALS — the generalisation of
/// [`recover_at_most_one`] to WIDER and MIXED-POLARITY counting cores (the ternary at-most-two cores of
/// parity-learning instances, and — via negation — at-least-`k`). Work in literal codes (`2·var + sign`,
/// negation flips the low bit). A width-`(k+1)` clause `ℓ₀ ∨ … ∨ ℓ_k` is violated only when all of
/// `{¬ℓ₀,…,¬ℓ_k}` are true, so it FORBIDS that literal set; a set `S` of literals (distinct variables) is an
/// at-most-`k` group iff EVERY `(k+1)`-subset of `S` is some clause's forbidden set. Greedy maximal extension
/// from each forbidden seed — a literal joins `S` only when every `k`-subset of `S` together with it is
/// forbidden, so the invariant keeps every `(k+1)`-subset forbidden and `at_most(S, k)` genuinely implied by
/// clauses present (sound to fuse). A clause `[¬a,¬b,¬c]` ⇒ at-most-2 of `{a,b,c}`; the three clauses
/// `[a,b],[a,c],[b,c]` ⇒ at-most-1 of `{¬a,¬b,¬c}` = **at-least-2** of `{a,b,c}`. `k = 1` over all-negative
/// binary clauses reproduces [`recover_at_most_one`].
pub fn recover_at_most_k(
    num_vars: usize,
    clauses: &[Vec<Lit_>],
    k: usize,
) -> Vec<crate::pseudo_boolean::PbConstraint> {
    use crate::pseudo_boolean::PbConstraint;
    use std::collections::HashSet;
    if k == 0 {
        return Vec::new();
    }
    let width = k + 1;
    let code = |l: &Lit_| (l.var() as usize) * 2 + usize::from(!l.is_positive()); // 2·var + sign
    let mut forbidden: HashSet<Vec<usize>> = HashSet::new();
    let mut cand_set: HashSet<usize> = HashSet::new();
    for c in clauses {
        if c.len() != width {
            continue;
        }
        let mut g: Vec<usize> = c.iter().map(|l| code(l) ^ 1).collect(); // the negated literals
        g.sort_unstable();
        let distinct_vars = g.windows(2).all(|w| w[0] >> 1 != w[1] >> 1);
        g.dedup();
        if g.len() == width && distinct_vars && g.iter().all(|&lc| lc >> 1 < num_vars) {
            for &lc in &g {
                cand_set.insert(lc);
            }
            forbidden.insert(g);
        }
    }
    if forbidden.is_empty() {
        return Vec::new();
    }
    let mut cand: Vec<usize> = cand_set.iter().copied().collect();
    cand.sort_unstable();
    let mut seeds: Vec<Vec<usize>> = forbidden.iter().cloned().collect();
    seeds.sort();
    let mut covered: HashSet<Vec<usize>> = HashSet::new();
    let mut out: Vec<PbConstraint> = Vec::new();
    // A combinatorial-work budget: each subset visit costs one. Stopping early is SOUND — a verified group
    // built so far (or any subset of it) is still a valid at-most-`k`, we just recover fewer / smaller ones.
    // It bounds the `C(group, k)` growth so a pathological wide-clause input can never blow up.
    let mut budget: usize = 1_000_000;
    const MAX_GROUP: usize = 96;
    for seed in seeds {
        if budget == 0 {
            break;
        }
        if covered.contains(&seed) {
            continue; // already inside a previously-emitted maximal group
        }
        let mut group = seed.clone();
        let mut group_vars: HashSet<usize> = group.iter().map(|&lc| lc >> 1).collect();
        for &v in &cand {
            if group_vars.contains(&(v >> 1)) || group.len() >= MAX_GROUP {
                continue; // distinct variables only (no `x` and `¬x` in one group)
            }
            // v joins iff every k-subset of `group` together with v is a forbidden (k+1)-subset.
            let mut starved = false;
            let ok = for_each_combo(&group, k, 0, &mut Vec::with_capacity(k), &mut |sub| {
                if budget == 0 {
                    starved = true;
                    return false;
                }
                budget -= 1;
                let mut s = sub.to_vec();
                s.push(v);
                s.sort_unstable();
                forbidden.contains(&s)
            });
            if starved {
                break;
            }
            if ok {
                group.push(v);
                group.sort_unstable();
                group_vars.insert(v >> 1);
            }
        }
        let mut fresh = false;
        for_each_combo(
            &group,
            width,
            0,
            &mut Vec::with_capacity(width),
            &mut |sub| {
                if budget == 0 {
                    return false;
                }
                budget -= 1;
                if covered.insert(sub.to_vec()) {
                    fresh = true;
                }
                true
            },
        );
        if fresh {
            let lits: Vec<(usize, bool)> = group.iter().map(|&lc| (lc >> 1, lc & 1 == 0)).collect();
            out.push(PbConstraint::at_most(&lits, k as i64));
        }
    }
    out
}

/// The widest core a fused solve recovers automatically: at-most-`k` groups up to this `k`. Beyond it,
/// call [`recover_at_most_k`] directly — but the `C(group, k)` scan cost (budget-bounded) and the
/// diminishing yield of ever-wider exclusion cliques make 4 (5-ary clauses) the default ceiling.
pub const MAX_RECOVERED_CARDINALITY: usize = 4;

/// Recover the cardinality substructure a fused solve consumes: every at-most-`k` group over arbitrary
/// literals, `k = 1 … `[`MAX_RECOVERED_CARDINALITY`] — pairwise, ternary, quaternary, … exclusion cliques of
/// any polarity, so at-most-`k`, at-least-`k` (the negated cliques), and mixed counting cores all flow into
/// the fused theory. Each [`recover_at_most_k`] returns immediately when the formula has no width-`(k+1)`
/// clause of the right shape, so widening is free on instances that lack those cores.
pub fn recover_cardinality_substructure(
    num_vars: usize,
    clauses: &[Vec<Lit_>],
) -> Vec<crate::pseudo_boolean::PbConstraint> {
    let mut out = Vec::new();
    for k in 1..=MAX_RECOVERED_CARDINALITY {
        out.extend(recover_at_most_k(num_vars, clauses, k));
    }
    out
}

/// **The fused parity + cardinality decision.** When a formula carries BOTH a GF(2) linear substructure
/// (recovered XOR equations via [`extract_xor`]) AND an at-most-one cardinality substructure (via
/// [`recover_at_most_one`]), decide it with the two live theories reasoning TOGETHER on one trail —
/// Gaussian elimination for the parity, cardinality propagation for the counting — the structural attack on
/// minimal-disagreement parity that neither alone cracks. Returns `Some(is_sat)`, or `None` when either
/// substructure is absent (so a caller falls through to other routes). Sound: the original clauses are
/// solved (Boolean-complete), the theories return only formula-entailed clauses, and a SAT model is
/// re-checked against the clauses (fail-closed). Uses the stateless [`crate::xor_engine::XorEngine`] (not
/// `IncXor`, whose value-blind trail-sync is unsafe without a backing clausal XOR encoding).
pub fn fused_parity_cardinality_decide(num_vars: usize, clauses: &[Vec<Lit_>]) -> Option<bool> {
    if num_vars == 0 {
        return None;
    }
    let eqs = extract_xor(num_vars, clauses);
    let amo = recover_cardinality_substructure(num_vars, clauses);
    if eqs.is_empty() || amo.is_empty() {
        return None;
    }
    // Break the ENTIRE symmetry group — permutations AND affine maps AND cross-compositions — COMPLETELY and
    // DYNAMICALLY via the aux-free SymmetryTheory, fused on the shared trail with parity + cardinality. No
    // static clauses, no aux variables.
    let mut s = Solver::new(num_vars);
    for c in clauses {
        s.add_clause(c.clone());
    }
    let mut theories: Vec<Box<dyn crate::cdcl::Theory>> = vec![
        Box::new(crate::xor_engine::XorEngine::new(num_vars, &eqs)),
        Box::new(crate::pseudo_boolean::CardinalityTheory::new(
            num_vars, &amo,
        )),
        Box::new(SymmetryTheory::new(
            num_vars,
            fused_symmetry_group(num_vars, clauses),
        )),
    ];
    match s.solve_with(&mut theories) {
        SolveResult::Sat(m) => clauses
            .iter()
            .all(|c| c.iter().any(|l| m[l.var() as usize] == l.is_positive()))
            .then_some(true),
        SolveResult::Unsat => Some(false),
    }
}

/// A **model-set (SEMANTIC) symmetry** checker for a fused instance: a variable permutation is a symmetry
/// iff it preserves the parity SOLUTION SPACE — the GF(2) row SPAN, not merely the gadget clauses, so it
/// sees the affine parity symmetry a syntactic clause-automorphism check is blind to (a permutation may map
/// one XOR equation to a *linear combination* of the others) — AND preserves the non-parity clauses as a
/// set. Sound: such a permutation maps models to models, so a lex-leader over it is satisfiability-preserving.
/// Built once per instance (the parity RREF basis is precomputed).
struct SemanticSymmetry {
    words: usize,
    rhs_bit: usize,
    basis: Vec<Vec<u64>>,
    pivots: Vec<usize>,
    eqs: Vec<(Vec<usize>, bool)>,
    non_parity: Vec<Vec<Lit_>>,
}

fn gf2_row(
    words: usize,
    rhs_bit: usize,
    vars: &[usize],
    rhs: bool,
    perm: Option<&[usize]>,
) -> Vec<u64> {
    let mut b = vec![0u64; words];
    for &v in vars {
        let idx = perm.map(|p| p[v]).unwrap_or(v);
        b[idx / 64] ^= 1 << (idx % 64);
    }
    if rhs {
        b[rhs_bit / 64] ^= 1 << (rhs_bit % 64);
    }
    b
}
fn gf2_reduce(v: &mut [u64], basis: &[Vec<u64>], pivots: &[usize]) {
    for (row, &piv) in basis.iter().zip(pivots) {
        if (v[piv / 64] >> (piv % 64)) & 1 == 1 {
            for w in 0..v.len() {
                v[w] ^= row[w];
            }
        }
    }
}
fn gf2_lowest(v: &[u64]) -> Option<usize> {
    v.iter()
        .enumerate()
        .find_map(|(w, &word)| (word != 0).then(|| w * 64 + word.trailing_zeros() as usize))
}

impl SemanticSymmetry {
    fn new(num_vars: usize, clauses: &[Vec<Lit_>]) -> Self {
        let eqs_raw = extract_xor(num_vars, clauses);
        let rhs_bit = num_vars;
        let words = num_vars.div_ceil(64) + 1;
        let eqs: Vec<(Vec<usize>, bool)> =
            eqs_raw.iter().map(|e| (e.vars.clone(), e.rhs)).collect();
        let mut xor_sets: std::collections::HashSet<Vec<usize>> = std::collections::HashSet::new();
        for (vars, _) in &eqs {
            let mut v = vars.clone();
            v.sort_unstable();
            v.dedup();
            xor_sets.insert(v);
        }
        let non_parity: Vec<Vec<Lit_>> = clauses
            .iter()
            .filter(|c| {
                let mut vs: Vec<usize> = c.iter().map(|l| l.var() as usize).collect();
                vs.sort_unstable();
                vs.dedup();
                !xor_sets.contains(&vs)
            })
            .cloned()
            .collect();
        let mut basis: Vec<Vec<u64>> = Vec::new();
        let mut pivots: Vec<usize> = Vec::new();
        for (vars, rhs) in &eqs {
            let mut v = gf2_row(words, rhs_bit, vars, *rhs, None);
            gf2_reduce(&mut v, &basis, &pivots);
            if let Some(p) = gf2_lowest(&v) {
                basis.push(v);
                pivots.push(p);
            }
        }
        SemanticSymmetry {
            words,
            rhs_bit,
            basis,
            pivots,
            eqs,
            non_parity,
        }
    }

    /// Is `perm` a model-set symmetry — parity span preserved AND non-parity clauses preserved?
    fn is_symmetry(&self, perm: &[usize]) -> bool {
        for (vars, rhs) in &self.eqs {
            let mut v = gf2_row(self.words, self.rhs_bit, vars, *rhs, Some(perm));
            gf2_reduce(&mut v, &self.basis, &self.pivots);
            if gf2_lowest(&v).is_some() {
                return false; // a permuted equation escapes the parity span
            }
        }
        let sigma = Perm::from_images(perm.iter().map(|&v| Lit_::pos(v as u32)).collect());
        perm_is_automorphism(&self.non_parity, &sigma)
    }
}

/// The **cardinality / parity symmetry seams.** Each recovered cardinality group is fully symmetric — its
/// members are interchangeable (`S_m` preserves "at most `k` of them true"). In a COUPLED instance only some
/// of those swaps survive as MODEL-SET symmetries; the parity coupling blocks the rest. This partitions the
/// candidate interchangeable pairs into `joint` (a transposition that preserves the model set — a free
/// symmetry to break) and `seams` (a transposition the parity genuinely tears). The check is SEMANTIC (parity
/// SPAN + non-parity clauses via `SemanticSymmetry`), so it recognizes the affine parity symmetry a
/// clause-automorphism test misses. The joint swaps generate the model-preserving symmetry at the seam.
pub struct CardinalitySeams {
    pub joint: Vec<(usize, usize)>,
    pub seams: Vec<(usize, usize)>,
}

/// Compute the [`CardinalitySeams`] of a fused instance: for every pair of variables interchangeable in some
/// recovered cardinality group, test whether swapping them is a model-set symmetry (`SemanticSymmetry`) —
/// joint if so, a seam if not. Pair-budget-bounded (a partial joint set still yields sound symmetry breaks).
pub fn cardinality_parity_seams(num_vars: usize, clauses: &[Vec<Lit_>]) -> CardinalitySeams {
    use std::collections::HashSet;
    const PAIR_BUDGET: usize = 20_000;
    let cons = recover_cardinality_substructure(num_vars, clauses);
    let checker = SemanticSymmetry::new(num_vars, clauses);
    let mut seen: HashSet<(usize, usize)> = HashSet::new();
    let mut pairs: Vec<(usize, usize)> = Vec::new();
    for pb in &cons {
        let mut vars: Vec<usize> = pb.terms().map(|(v, _, _)| v).collect();
        vars.sort_unstable();
        vars.dedup();
        for i in 0..vars.len() {
            for j in (i + 1)..vars.len() {
                if seen.insert((vars[i], vars[j])) {
                    pairs.push((vars[i], vars[j]));
                }
            }
        }
    }
    pairs.sort_unstable();
    let mut joint = Vec::new();
    let mut seams = Vec::new();
    for &(a, b) in pairs.iter().take(PAIR_BUDGET) {
        let mut perm: Vec<usize> = (0..num_vars).collect();
        perm.swap(a, b);
        if checker.is_symmetry(&perm) {
            joint.push((a, b));
        } else {
            seams.push((a, b));
        }
    }
    CardinalitySeams { joint, seams }
}

/// Break the joint cardinality/parity symmetry: union the joint swaps into interchangeable components (a
/// component is fully symmetric — the group generated by its transpositions is `S`, every element an
/// automorphism), then order each component's variables into a descending chain with `vᵢ ≥ vᵢ₊₁` clauses
/// `(vᵢ ∨ ¬vᵢ₊₁)`. A lex-leader: equisatisfiable (every orbit keeps its ordered representative), so it
/// shrinks the fused search without changing the verdict. Returns the extra clauses (empty if no joint
/// symmetry).
pub fn cardinality_symmetry_break(num_vars: usize, clauses: &[Vec<Lit_>]) -> Vec<Vec<Lit_>> {
    use std::collections::BTreeSet;
    let joint = cardinality_parity_seams(num_vars, clauses).joint;
    if joint.is_empty() {
        return Vec::new();
    }
    let mut parent: Vec<usize> = (0..num_vars).collect();
    fn find(parent: &mut [usize], mut x: usize) -> usize {
        while parent[x] != x {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        x
    }
    for &(a, b) in &joint {
        let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
        if ra != rb {
            parent[ra] = rb;
        }
    }
    let mut members: std::collections::HashMap<usize, BTreeSet<usize>> =
        std::collections::HashMap::new();
    for &(a, b) in &joint {
        let r = find(&mut parent, a);
        let set = members.entry(r).or_default();
        set.insert(a);
        set.insert(b);
    }
    let mut comps: Vec<Vec<usize>> = members
        .into_values()
        .map(|s| s.into_iter().collect())
        .collect();
    comps.sort();
    let mut out = Vec::new();
    for vs in comps {
        for w in vs.windows(2) {
            out.push(vec![Lit_::pos(w[0] as u32), Lit_::neg(w[1] as u32)]); // v_i ≥ v_{i+1}
        }
    }
    out
}

/// The generators of the joint parity/cardinality symmetry, UP THE WREATH CHAIN. Two levels:
/// * **within-orbit** — the joint swaps (interchangeable variables), the *class* `S_m` inside each orbit;
/// * **block** — whole joint orbits mapped to each other across the seams (`(0 2)(1 3)` when the individual
///   `0↔2`, `1↔3` are seams): automorphisms of the quotient, the *family* `S_k` permuting the `k` blocks.
///
/// Together these generate the wreath product `S_m ≀ S_k`. Feed to [`crate::sym_break::lex_leader_sbp`] to
/// break class and family at once — pulling the symmetry apart on the seams, not just within them. Block
/// pairs are tested under the sorted correspondence (a sound heuristic — only genuine automorphisms are
/// emitted), budget-bounded.
pub fn cardinality_symmetry_generators(num_vars: usize, clauses: &[Vec<Lit_>]) -> Vec<Vec<usize>> {
    let seams = cardinality_parity_seams(num_vars, clauses);
    // A variable permutation is a MODEL-SET symmetry (parity span + non-parity clauses) — the semantic check,
    // so the block/wreath levels recognize the affine parity symmetry a clause-automorphism test misses.
    let checker = SemanticSymmetry::new(num_vars, clauses);
    let is_auto = |perm: &[usize]| -> bool { checker.is_symmetry(perm) };
    // within-orbit swaps — already verified joint, so pushed directly as transposition permutations.
    let mut gens: Vec<Vec<usize>> = Vec::new();
    for &(a, b) in &seams.joint {
        let mut p: Vec<usize> = (0..num_vars).collect();
        p.swap(a, b);
        gens.push(p);
    }

    // Joint orbits (union-find over the within-orbit swaps).
    let mut parent: Vec<usize> = (0..num_vars).collect();
    fn find(p: &mut [usize], mut x: usize) -> usize {
        while p[x] != x {
            p[x] = p[p[x]];
            x = p[x];
        }
        x
    }
    for &(a, b) in &seams.joint {
        let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
        if ra != rb {
            parent[ra] = rb;
        }
    }
    let mut orbmap: std::collections::HashMap<usize, std::collections::BTreeSet<usize>> =
        std::collections::HashMap::new();
    for &(a, b) in &seams.joint {
        let r = find(&mut parent, a);
        let set = orbmap.entry(r).or_default();
        set.insert(a);
        set.insert(b);
    }
    let mut orbits: Vec<Vec<usize>> = orbmap
        .into_values()
        .map(|s| s.into_iter().collect())
        .collect();
    orbits.sort();

    // Climb the wreath chain to the TOP: each round, find interchangeable blocks (same length, position-wise
    // swap an automorphism), emit a spanning set of block swaps, then MERGE each interchange-class into one
    // larger block (its members concatenated in canonical order) so the next round sees the symmetry of the
    // QUOTIENT. Repeat until a round finds no interchangeable pair (the top). Budget-bounded; sound —
    // merging is exact (a super-block swap maps sub-block k to sub-block k because both keep canonical order).
    let mut blocks: Vec<Vec<usize>> = orbits;
    let mut budget: usize = 4000;
    loop {
        let n = blocks.len();
        if n < 2 || budget == 0 {
            break;
        }
        let mut bp: Vec<usize> = (0..n).collect();
        let mut level_gens: Vec<Vec<usize>> = Vec::new();
        'lvl: for i in 0..n {
            for j in (i + 1)..n {
                if budget == 0 {
                    break 'lvl;
                }
                if blocks[i].len() != blocks[j].len() {
                    continue;
                }
                let (ri, rj) = (find(&mut bp, i), find(&mut bp, j));
                if ri == rj {
                    continue; // already merged this round
                }
                budget -= 1;
                let mut p: Vec<usize> = (0..num_vars).collect();
                for (&a, &b) in blocks[i].iter().zip(&blocks[j]) {
                    p.swap(a, b);
                }
                if is_auto(&p) {
                    bp[ri] = rj; // spanning-tree edge
                    level_gens.push(p);
                }
            }
        }
        if level_gens.is_empty() {
            break; // top of the chain
        }
        gens.extend(level_gens);
        let mut classes: std::collections::BTreeMap<usize, Vec<usize>> =
            std::collections::BTreeMap::new();
        for i in 0..n {
            classes.entry(find(&mut bp, i)).or_default().push(i);
        }
        let mut next: Vec<Vec<usize>> = Vec::new();
        for (_, idxs) in classes {
            let mut mem: Vec<Vec<usize>> = idxs.iter().map(|&i| blocks[i].clone()).collect();
            mem.sort_by_key(|b| b[0]);
            next.push(mem.into_iter().flatten().collect());
        }
        next.sort();
        blocks = next;
    }
    gens
}

/// Detect **affine (shear) symmetries** of a fused instance that preserve the model set — the affine parity
/// symmetries a variable-permutation break structurally cannot express (an image bit becomes an XOR of two
/// variables). For a variable `i` that appears in NO non-parity (cardinality / residual) clause and any
/// `j ≠ i`, the shear `x_i ↦ x_i ⊕ x_j` moves only position `i`, so it leaves the cardinality/residual
/// intact; it is a model-set symmetry iff it maps the parity SOLUTION SPACE to itself (checked on the
/// space's affine spanning set — particular ⊕ each kernel basis vector). Each valid shear is returned as an
/// SBP affine-map spec (identity except output `i = ⊕{i,j}`) for [`crate::sym_break::affine_lex_leader_sbp`].
/// Gated to `num_vars ≤ 64`.
pub fn affine_parity_symmetries(
    num_vars: usize,
    clauses: &[Vec<Lit_>],
) -> Vec<Vec<(Vec<usize>, bool)>> {
    if num_vars == 0 || num_vars > 64 {
        return Vec::new();
    }
    let eqs = extract_xor(num_vars, clauses);
    if eqs.is_empty() {
        return Vec::new();
    }
    let rows: Vec<u64> = eqs
        .iter()
        .map(|e| e.vars.iter().fold(0u64, |a, &v| a | (1u64 << v)))
        .collect();
    let rhs: Vec<bool> = eqs.iter().map(|e| e.rhs).collect();
    let Some(space) = crate::gf2::solve_gf2(num_vars, &rows, &rhs) else {
        return Vec::new(); // inconsistent parity — no solution space
    };
    let mut xor_sets: std::collections::HashSet<Vec<usize>> = std::collections::HashSet::new();
    for e in &eqs {
        let mut v = e.vars.clone();
        v.sort_unstable();
        v.dedup();
        xor_sets.insert(v);
    }
    let mut moved_forbidden = vec![false; num_vars]; // vars a shear must not move (in a non-parity clause)
    for c in clauses {
        let mut vs: Vec<usize> = c.iter().map(|l| l.var() as usize).collect();
        vs.sort_unstable();
        vs.dedup();
        if !xor_sets.contains(&vs) {
            for &v in &vs {
                moved_forbidden[v] = true;
            }
        }
    }
    let mut points: Vec<Vec<bool>> = vec![space.particular.clone()];
    for k in &space.kernel_basis {
        let mut p = space.particular.clone();
        for v in 0..num_vars {
            p[v] ^= k[v];
        }
        points.push(p);
    }
    let satisfies = |x: &[bool]| {
        eqs.iter()
            .all(|e| e.vars.iter().fold(false, |a, &v| a ^ x[v]) == e.rhs)
    };
    // Build (and VERIFY, sound regardless of the algebra) one affine generator: `moved` = the coordinates it
    // changes (must all be parity-only); `source = Some(j)` is the transvection `x ↦ x ⊕ x_j·1_moved`,
    // `source = None` is the translation `x ↦ x ⊕ 1_moved`. Kept only if it maps the solution space to itself.
    let make_map = |moved: &[usize], source: Option<usize>| -> Option<Vec<(Vec<usize>, bool)>> {
        if moved.is_empty() || moved.iter().any(|&i| moved_forbidden[i]) {
            return None;
        }
        let preserves = points.iter().all(|x| {
            let add = source.map_or(true, |j| x[j]);
            if !add {
                return true;
            }
            let mut y = x.clone();
            for &i in moved {
                y[i] ^= true;
            }
            satisfies(&y)
        });
        if !preserves {
            return None;
        }
        let mut spec: Vec<(Vec<usize>, bool)> = (0..num_vars).map(|k| (vec![k], false)).collect();
        for &i in moved {
            spec[i] = match source {
                Some(j) => (vec![i, j], false),
                None => (vec![i], true),
            };
        }
        Some(spec)
    };
    let mut maps: Vec<Vec<(Vec<usize>, bool)>> = Vec::new();
    // Single-variable shears x_i ↦ x_i ⊕ x_j.
    for i in 0..num_vars {
        for j in 0..num_vars {
            if j != i {
                if let Some(m) = make_map(&[i], Some(j)) {
                    maps.push(m);
                }
            }
        }
    }
    // The GL rung: multi-coordinate transvections and translations along `K ∩ P` — the kernel directions of
    // the parity supported entirely on parity-only variables. For each such direction κ: the translation
    // `x ↦ x ⊕ κ` (flip all of κ's support together) and the transvections `x ↦ x ⊕ x_j·κ` (add x_j to all of
    // κ's support). These mix SEVERAL coupled parity variables at once — the affine symmetries neither a
    // single shear nor any permutation expresses.
    for kappa in kernel_intersect_p(&space.kernel_basis, &moved_forbidden, num_vars) {
        let support: Vec<usize> = (0..num_vars).filter(|&i| kappa[i]).collect();
        if let Some(m) = make_map(&support, None) {
            maps.push(m);
        }
        for j in 0..num_vars {
            if !kappa[j] {
                if let Some(m) = make_map(&support, Some(j)) {
                    maps.push(m);
                }
            }
        }
    }
    maps.sort();
    maps.dedup();
    maps
}

/// `K ∩ P`: the parity kernel directions supported entirely on parity-only variables (zero on every
/// `moved_forbidden` coordinate). Solve, over the kernel-basis coefficients `a`, the homogeneous system "the
/// combination `Σ aₜ κₜ` is zero on each forbidden coordinate" ([`crate::gf2::solve_gf2`]); each null-space
/// vector `a` yields the direction `Σ aₜ κₜ`. Empty forbidden set ⇒ all of `K`.
fn kernel_intersect_p(
    kernel_basis: &[Vec<bool>],
    moved_forbidden: &[bool],
    num_vars: usize,
) -> Vec<Vec<bool>> {
    let d = kernel_basis.len();
    if d == 0 || d > 64 {
        return Vec::new();
    }
    let forbidden: Vec<usize> = (0..num_vars).filter(|&c| moved_forbidden[c]).collect();
    let rows: Vec<u64> = forbidden
        .iter()
        .map(|&c| {
            (0..d).fold(0u64, |acc, t| {
                if kernel_basis[t][c] {
                    acc | (1u64 << t)
                } else {
                    acc
                }
            })
        })
        .collect();
    let rhs = vec![false; rows.len()];
    let Some(null) = crate::gf2::solve_gf2(d, &rows, &rhs) else {
        return Vec::new();
    };
    null.kernel_basis
        .iter()
        .map(|a| {
            let mut kappa = vec![false; num_vars];
            for (t, &at) in a.iter().enumerate().take(d) {
                if at {
                    for v in 0..num_vars {
                        kappa[v] ^= kernel_basis[t][v];
                    }
                }
            }
            kappa
        })
        .collect()
}

/// **Every** single-transposition model-set symmetry — every variable swap that preserves the model set
/// (`SemanticSymmetry`): parity-variable permutations, cardinality permutations, and residual-clause
/// permutations alike. Broadens the wreath-tower (cardinality-only) symmetry to the FULL permutation
/// symmetry of the whole instance — in particular the permutations of parity-only variables neither the
/// cardinality detector nor the affine shear detector finds. Budget-bounded (`n ≤ 64`, `C(n,2)` checks).
pub fn all_transposition_symmetries(num_vars: usize, clauses: &[Vec<Lit_>]) -> Vec<Vec<usize>> {
    if num_vars < 2 || num_vars > 64 {
        return Vec::new();
    }
    let checker = SemanticSymmetry::new(num_vars, clauses);
    let mut gens: Vec<Vec<usize>> = Vec::new();
    const BUDGET: usize = 4096;
    let mut checks = 0usize;
    'a: for a in 0..num_vars {
        for b in (a + 1)..num_vars {
            if checks >= BUDGET {
                break 'a;
            }
            checks += 1;
            let mut p: Vec<usize> = (0..num_vars).collect();
            p.swap(a, b);
            if checker.is_symmetry(&p) {
                gens.push(p);
            }
        }
    }
    gens
}

/// An **affine map over GF(2)** on `n ≤ 64` variables: `α(x)[j] = parity(rows[j] & x) ⊕ ((trans >> j) & 1)`.
/// Permutations, shears, translations — every symmetry we detect — are affine, so unifying them here lets us
/// compose them into ONE group and break it completely.
#[derive(Clone, PartialEq, Eq, Hash)]
struct AffineMap {
    rows: Vec<u64>,
    trans: u64,
}
impl AffineMap {
    fn identity(n: usize) -> Self {
        AffineMap {
            rows: (0..n).map(|j| 1u64 << j).collect(),
            trans: 0,
        }
    }
    fn from_perm(p: &[usize]) -> Self {
        AffineMap {
            rows: p.iter().map(|&pj| 1u64 << pj).collect(),
            trans: 0,
        }
    }
    fn from_spec(spec: &[(Vec<usize>, bool)], n: usize) -> Self {
        let mut rows = vec![0u64; n];
        let mut trans = 0u64;
        for (j, (xs, b)) in spec.iter().enumerate().take(n) {
            rows[j] = xs.iter().fold(0u64, |a, &v| a | (1u64 << v));
            if *b {
                trans |= 1u64 << j;
            }
        }
        AffineMap { rows, trans }
    }
    fn is_identity(&self) -> bool {
        self.trans == 0 && self.rows.iter().enumerate().all(|(j, &r)| r == (1u64 << j))
    }
    /// `self ∘ other` (apply `other` first). `(A∘B).rows[j] = ⊕_{k∈A.rows[j]} B.rows[k]`,
    /// `(A∘B).trans[j] = A.trans[j] ⊕ ⊕_{k∈A.rows[j]} B.trans[k]`.
    fn compose(&self, other: &AffineMap) -> AffineMap {
        let n = self.rows.len();
        let mut rows = vec![0u64; n];
        let mut trans = 0u64;
        for j in 0..n {
            let (mut r, mut t) = (0u64, (self.trans >> j) & 1);
            let mut sr = self.rows[j];
            while sr != 0 {
                let k = sr.trailing_zeros() as usize;
                sr &= sr - 1;
                r ^= other.rows[k];
                t ^= (other.trans >> k) & 1;
            }
            rows[j] = r;
            trans |= t << j;
        }
        AffineMap { rows, trans }
    }
    fn to_spec(&self) -> Vec<(Vec<usize>, bool)> {
        self.rows
            .iter()
            .enumerate()
            .map(|(j, &r)| {
                (
                    (0..64).filter(|&b| (r >> b) & 1 == 1).collect(),
                    (self.trans >> j) & 1 == 1,
                )
            })
            .collect()
    }
}

/// Close a set of affine generators into their group by BFS (right multiplication). Our generators are
/// involutions (transpositions, shears, translations), so right-multiplication reaches the whole group.
/// Returns `None` if the closure exceeds `cap` (the caller falls back to a partial, generator-only break).
fn affine_group_closure(gens: &[AffineMap], num_vars: usize, cap: usize) -> Option<Vec<AffineMap>> {
    use std::collections::HashSet;
    let id = AffineMap::identity(num_vars);
    let mut seen: HashSet<AffineMap> = HashSet::from([id.clone()]);
    let mut frontier = vec![id];
    while let Some(g) = frontier.pop() {
        for gen in gens {
            let h = g.compose(gen);
            if seen.insert(h.clone()) {
                if seen.len() > cap {
                    return None;
                }
                frontier.push(h);
            }
        }
    }
    Some(seen.into_iter().collect())
}

/// **The ultimate symmetry break for the fused route** — the ENTIRE symmetry group, broken COMPLETELY and
/// DYNAMICALLY with ZERO aux variables. Every symmetry we detect — the wreath-tower permutations (class ×
/// family), every transposition symmetry, AND the affine parity maps (shears / GL-part transvections /
/// translations) with their cross-compositions — is affine over GF(2), so we unify them into ONE generating
/// set, close it into the FULL GROUP, and hand every element to the aux-free [`SymmetryTheory`], which
/// enforces `x ≤_lex α(x)` for each `α` by propagation on the shared trail. That is the COMPLETE break
/// (exactly one representative per orbit) of the full permutation × affine group, with no static clauses and
/// no variable blow-up. When the group exceeds the cap (or `num_vars > 64`) it degrades to the generators
/// (partial, still sound — ≥ 1 representative per orbit). Returns the group as affine-map specs for the theory.
/// Convert a signed-permutation clause automorphism `σ` (a [`crate::proof::Perm`], `σ(+v)` a literal) into
/// an affine spec. The assignment action is `y(σ(l)) = x(l)`, so for `σ(+v) = ±w` the output bit `w` reads
/// source `v` (with a translation bit when the image is negative): `σ(+v)=+w ⇒ y_w = x_v`,
/// `σ(+v)=-w ⇒ y_w = ¬x_v`. Initialized to the identity so a fixed variable stays fixed.
fn signed_perm_to_spec(p: &crate::proof::Perm, num_vars: usize) -> Vec<(Vec<usize>, bool)> {
    let mut spec: Vec<(Vec<usize>, bool)> = (0..num_vars).map(|w| (vec![w], false)).collect();
    for v in 0..num_vars {
        let img = p.apply(Lit::pos(v as u32));
        let w = img.var() as usize;
        if w < num_vars {
            spec[w] = (vec![v], !img.is_positive());
        }
    }
    spec
}

pub fn fused_symmetry_group(
    num_vars: usize,
    clauses: &[Vec<Lit_>],
) -> Vec<Vec<(Vec<usize>, bool)>> {
    let perm_gens = fused_permutation_generators(num_vars, clauses);
    if num_vars < 1 || num_vars > 64 {
        // The affine machinery (u64 rows) is off — the permutation generators as specs (partial break).
        return perm_gens
            .iter()
            .map(|p| p.iter().map(|&pi| (vec![pi], false)).collect())
            .collect();
    }
    let aff_specs = affine_parity_symmetries(num_vars, clauses);
    // The FULL syntactic clause-automorphism generators — individualization-refinement (bliss/saucy-style),
    // catching multi-cycle and SIGNED (negation) generators that the transposition/affine detectors miss.
    // A clause automorphism is a model-set symmetry, so its lex-leader break is equisatisfiable.
    let syntactic = crate::symmetry_detect::find_generators(num_vars, clauses);
    let mut gens: Vec<AffineMap> = perm_gens.iter().map(|p| AffineMap::from_perm(p)).collect();
    gens.extend(aff_specs.iter().map(|s| AffineMap::from_spec(s, num_vars)));
    gens.extend(
        syntactic
            .iter()
            .map(|p| AffineMap::from_spec(&signed_perm_to_spec(p, num_vars), num_vars)),
    );
    gens.retain(|g| !g.is_identity());
    if gens.is_empty() {
        return Vec::new();
    }
    const CAP: usize = 2048;
    match affine_group_closure(&gens, num_vars, CAP) {
        Some(group) => group
            .iter()
            .filter(|g| !g.is_identity())
            .map(|g| g.to_spec())
            .collect(),
        None => gens.iter().map(|g| g.to_spec()).collect(), // group too large — the generators (partial)
    }
}

/// The permutation generators for the DYNAMIC symmetry theory: the wreath-tower cardinality permutations
/// plus every transposition symmetry of the whole instance.
pub fn fused_permutation_generators(num_vars: usize, clauses: &[Vec<Lit_>]) -> Vec<Vec<usize>> {
    let mut g = cardinality_symmetry_generators(num_vars, clauses);
    g.extend(all_transposition_symmetries(num_vars, clauses));
    g
}

/// A live **symmetry theory** for [`crate::cdcl::Solver::solve_with`] — DYNAMIC lex-leader propagation over
/// the whole symmetry group. Each element is an AFFINE map `α` (a model-set automorphism) given as a
/// per-output spec `α(x)[i] = ⊕_{s∈xset_i} x_s ⊕ b_i`; permutations are the special case where every `xset`
/// is a single variable. During search it enforces `x ≤_lex α(x)`: walking positions while the prefix stays
/// equal, computing `α(x)[i]` from the trail (an XOR of the support), then at the frontier forcing `x_i = 0`
/// when `α(x)[i] = 0`, forcing the last free support bit to make `α(x)[i] = 1` when `x_i = 1`, or conflicting
/// — AUX-FREE, on the shared trail, fused alongside the parity and cardinality theories. Sound: the orbit's
/// global lex-minimum satisfies `x ≤_lex α(x)` for EVERY `α`, so the whole group's predicate keeps exactly
/// one representative per orbit (equisatisfiable), and every reason clause (prefix + support witnesses) is
/// implied by it.
pub struct SymmetryTheory {
    num_vars: usize,
    maps: Vec<Vec<(Vec<usize>, bool)>>,
}
impl SymmetryTheory {
    pub fn new(num_vars: usize, maps: Vec<Vec<(Vec<usize>, bool)>>) -> Self {
        SymmetryTheory { num_vars, maps }
    }
    /// Convenience: build from variable permutations (each `σ[i]` the image variable of output `i`).
    pub fn from_perms(num_vars: usize, perms: Vec<Vec<usize>>) -> Self {
        let maps = perms
            .into_iter()
            .map(|p| p.into_iter().map(|pi| (vec![pi], false)).collect())
            .collect();
        SymmetryTheory { num_vars, maps }
    }
}
impl crate::cdcl::Theory for SymmetryTheory {
    fn propagate(&mut self, trail: &[crate::cdcl::Lit]) -> Vec<Vec<crate::cdcl::Lit>> {
        use crate::cdcl::Lit;
        let n = self.num_vars;
        let mut a: Vec<Option<bool>> = vec![None; n];
        for &l in trail {
            let v = l.var() as usize;
            if v < n {
                a[v] = Some(l.is_positive());
            }
        }
        // The assigned support literals of an output, each FALSE now — they witness the value it evaluated to.
        let support_witness = |xset: &[usize], a: &[Option<bool>]| -> Vec<Lit> {
            xset.iter()
                .filter_map(|&s| a[s].map(|sv| Lit::new(s as u32, !sv)))
                .collect()
        };
        let mut out = Vec::new();
        for map in &self.maps {
            let mut prefix: Vec<Lit> = Vec::new();
            for i in 0..n {
                let (xset, b) = &map[i];
                if xset.len() == 1 && xset[0] == i && !*b {
                    continue; // identity output — α(x)[i] = x_i, always equal
                }
                // Evaluate α(x)[i] = ⊕ support ⊕ b, tracking the unassigned support.
                let mut val = *b;
                let mut free: Vec<usize> = Vec::new();
                for &s in xset {
                    match a[s] {
                        Some(sv) => val ^= sv,
                        None => free.push(s),
                    }
                }
                match (a[i], free.len()) {
                    (Some(vi), 0) if vi == val => {
                        // x_i = α(x)[i]: prefix stays equal — witness both.
                        prefix.push(Lit::new(i as u32, !vi));
                        prefix.extend(support_witness(xset, &a));
                    }
                    (Some(false), 0) => break, // x_i=0 < α(x)[i]=1 ⇒ lex satisfied
                    // α(x)[i]=0 (val is false here, since the equal case was handled above): conflict when
                    // x_i=1, force x_i=0 when free. Clause {prefix} ∨ {¬support} ∨ ¬x_i.
                    (Some(true), 0) => {
                        let mut c = prefix.clone();
                        c.extend(support_witness(xset, &a));
                        c.push(Lit::new(i as u32, false));
                        out.push(c);
                        break;
                    }
                    (None, 0) if !val => {
                        let mut c = prefix.clone();
                        c.extend(support_witness(xset, &a));
                        c.push(Lit::new(i as u32, false));
                        out.push(c);
                        break;
                    }
                    (Some(true), 1) => {
                        // x_i=1, one free support bit s*: force it so α(x)[i]=1. `val` is the XOR of the
                        // assigned support ⊕ b, so x_{s*} must be ¬val.
                        let s = free[0];
                        let mut c = prefix.clone();
                        c.push(Lit::new(i as u32, false)); // ¬x_i (false now)
                        for &o in xset {
                            if o != s {
                                c.push(Lit::new(o as u32, !a[o].unwrap()));
                            }
                        }
                        c.push(Lit::new(s as u32, !val)); // force x_{s*} = ¬val
                        out.push(c);
                        break;
                    }
                    _ => break, // frontier undetermined
                }
            }
        }
        out
    }
}

/// The result of the unified engine: which physics collapsed the formula, with its checkable
/// artifact.
#[derive(Clone, Debug)]
pub enum AutoCollapse {
    /// Geometric collapse — covering symmetry, with the discovered measure and certified refutation.
    Geometric {
        measure: CollapsingMeasure,
        ranked: RankedRefutation,
    },
    /// Cardinality collapse — a covering with more items than bins refuted by cutting planes (no
    /// symmetry required), with the descent trajectory and the constraint count.
    Cardinality {
        trajectory: Vec<u64>,
        reached_goal: bool,
        constraints: usize,
    },
    /// Algebraic collapse — parity, with the GF(2) Lyapunov trajectory and the XOR count.
    Algebraic {
        trajectory: Vec<u64>,
        reached_goal: bool,
        xor_equations: usize,
    },
    /// No collapse found in any class — a bounded impossibility (the formula may still be UNSAT;
    /// no covering, cardinality, or parity structure was recognized).
    None,
}

/// **The unified, structure-recognizing engine.** Given only opaque clauses, decide WHICH collapse
/// physics applies — covering symmetry or parity — and produce the corresponding Lyapunov-certified
/// artifact. One engine; the structure selects the mechanism. Sound and fail-closed on both paths.
pub fn auto_collapse(num_vars: usize, formula: &[Vec<Lit_>]) -> AutoCollapse {
    // 1. Geometric: is there a covering *symmetry*? Preferred — it yields a checkable PR/SR proof.
    //    (fail-closed — `ranked.refuted` only if checked).
    if let Some((measure, ranked)) = solve_by_measure_synthesis(num_vars, formula) {
        if ranked.refuted {
            return AutoCollapse::Geometric { measure, ranked };
        }
    }
    // 2. Cardinality: a covering with no usable symmetry can still collapse by cutting planes
    //    (more items than bins ⇒ 0 ≥ 1). This catches *asymmetric* coverings the geometric route
    //    misses — e.g. the mutilated chessboard.
    if let Some((trajectory, reached_goal, constraints)) = cardinality_collapse(num_vars, formula) {
        if reached_goal {
            return AutoCollapse::Cardinality {
                trajectory,
                reached_goal,
                constraints,
            };
        }
    }
    // 3. Algebraic: recover the parity structure and collapse it over GF(2).
    let eqs = extract_xor(num_vars, formula);
    if !eqs.is_empty() {
        let (trajectory, reached_goal) = gaussian_lyapunov(&eqs, num_vars);
        if reached_goal {
            return AutoCollapse::Algebraic {
                trajectory,
                reached_goal,
                xor_equations: eqs.len(),
            };
        }
    }
    AutoCollapse::None
}

#[cfg(test)]
#[path = "../../logicaffeine_tests/unit/logicaffeine_proof/lyapunov/tests.rs"]
mod tests;