symplex 0.9.0

//! Exact linear programming over ℚ: two-phase simplex with Bland's rule,
//! dual values and Farkas infeasibility certificates.
//!
//! Everything is exact — no tolerances, no "numerically infeasible"
//! verdicts.  Problem data, optima, duals and certificates are
//! `Ratio<BigInt>`; internally the dense tableau pivots on integers with a
//! common denominator (see [Algorithm and cost](#algorithm-and-cost)).
//! Use [`linprog_matrix`] to feed [`Matrix`] data whose entries are
//! numeric literals.
//!
//! # Problem form
//!
//! [`LpProblem`] is a builder for
//!
//! ```text
//! minimise / maximise   cᵀx
//! subject to            aᵢ·x  ≤ / = / ≥  bᵢ        (one row per constraint)
//!                       lⱼ ≤ xⱼ ≤ uⱼ                 (default 0 ≤ xⱼ < ∞)
//! ```
//!
//! [`linprog`] is a SciPy-shaped convenience wrapper; [`feasible_nonneg`]
//! answers "is there an `x ≥ 0` with `A x = b`?" exactly,
//! [`feasible_nonneg_certified`] additionally returns the Farkas certificate
//! when the answer is no, and [`nonneg_combination`] asks the same question
//! about a target vector and a list of generating vectors (cone membership).
//!
//! # Algorithm and cost
//!
//! Two-phase dense simplex.  Pivots follow Dantzig's rule; after twelve
//! consecutive degenerate (zero-length) steps Bland's rule takes over until
//! the next improving step, so cycling is impossible while the highly
//! degenerate certificate LPs are not condemned to Bland's slow walk from
//! their first zero pivot.  Every `≤`/`≥` row gets a
//! slack, every row gets an artificial (their columns double as `B⁻¹`, from
//! which the duals are read), free variables are split `x = x⁺ − x⁻`, finite
//! lower bounds are shifted away and a finite *upper* bound on a variable
//! that also has a lower bound costs one extra row.
//!
//! The tableau is kept **integral** (Bareiss/Edmonds integer pivoting):
//! each constraint row is scaled once by the lcm of its denominators, and
//! every pivot applies the fraction-free rule `row ← (row·p − row[s]·pivot_row)
//! / d`, so all entries are integers sharing the common denominator `d`
//! (the current pivot, `±det B`).  Entering and leaving variables are
//! chosen from the same rational values a `Ratio` tableau would hold — the
//! pivot sequence, optimum, duals and certificates are identical — but no
//! gcd is computed in the pivot loop.  Each pivot is `O(m·n)` big-integer
//! operations; a 60-row × 160-variable problem with default bounds solves
//! in under 100 ms in a release build.
//!
//! # Duals and certificates
//!
//! When the status is [`LpStatus::Optimal`], [`LpSolution::duals`] holds one
//! **shadow price** `yᵢ` per constraint, in insertion order: the rate of
//! change of the optimal objective value *of the problem as posed* with
//! respect to `bᵢ`.  Consequently, for a minimisation `yᵢ ≤ 0` on `≤` rows
//! and `yᵢ ≥ 0` on `≥` rows (the signs flip for a maximisation), `yᵢ` is
//! free on `=` rows, and with the reduced costs `r = c − Aᵀy`:
//!
//! * complementary slackness: `yᵢ·(aᵢ·x* − bᵢ) = 0` for every row;
//! * `rⱼ = 0` unless `x*ⱼ` sits at a finite bound (for a minimisation
//!   `rⱼ ≥ 0` at a lower bound and `rⱼ ≤ 0` at an upper bound; reversed for
//!   a maximisation);
//! * strong duality: `cᵀx* = yᵀb + Σⱼ rⱼ x*ⱼ`, which reduces to
//!   `cᵀx* = yᵀb` under the default bounds `x ≥ 0`.
//!
//! When the status is [`LpStatus::Infeasible`], [`LpSolution::farkas`] is a
//! **Farkas certificate** `y` (one entry per constraint) with `yᵢ ≥ 0` on
//! `≤` rows, `yᵢ ≤ 0` on `≥` rows, free on `=` rows, such that, with
//! `g = Aᵀy`,
//!
//! ```text
//! Σⱼ  inf { gⱼ·xⱼ : lⱼ ≤ xⱼ ≤ uⱼ }   >   yᵀb
//! ```
//!
//! where every infimum is finite (`gⱼ > 0 ⇒ lⱼ` finite, `gⱼ < 0 ⇒ uⱼ`
//! finite, `gⱼ = 0` contributes `0`).  Any feasible `x` would satisfy
//! `(Aᵀy)·x ≤ yᵀb`, so the inequality proves that none exists.  With no
//! finite bounds this is the textbook form `Aᵀy = 0, yᵀb < 0`.  The
//! certificate is `None` only when infeasibility is caused by the bounds
//! alone (`lⱼ > uⱼ`).
//!
//! # Examples
//!
//! ```
//! use symplex::linprog::{LpProblem, LpStatus, qi};
//!
//! // max 3x + 2y  s.t.  x + y ≤ 4,  x + 3y ≤ 6,  x, y ≥ 0
//! let sol = LpProblem::maximize(vec![qi(3), qi(2)])
//!     .le(vec![qi(1), qi(1)], qi(4))
//!     .le(vec![qi(1), qi(3)], qi(6))
//!     .solve()
//!     .unwrap();
//! assert_eq!(sol.status, LpStatus::Optimal);
//! assert_eq!(sol.x, vec![qi(4), qi(0)]);
//! assert_eq!(sol.objective, Some(qi(12)));
//! // Shadow prices: the first constraint is binding with price 3.
//! assert_eq!(sol.duals, vec![qi(3), qi(0)]);
//! ```

use std::fmt;

use num_bigint::BigInt;
use num_integer::Integer;
use num_rational::Ratio;
use num_traits::{One, Signed, Zero};

use crate::api::context::Context;
use crate::api::expr::Ex;
use crate::base::errors::SymplexError;
use crate::domains::matrix::Matrix;

/// Exact rational number used throughout this module.
pub type Q = Ratio<BigInt>;

/// The rational `n / d`.
///
/// # Panics
///
/// Panics if `d == 0` (a programming error, like a zero literal
/// denominator).
///
/// # Examples
///
/// ```
/// use symplex::linprog::q;
/// assert_eq!(q(2, 4), q(1, 2));
/// ```
pub fn q(n: i64, d: i64) -> Q {
    Ratio::new(BigInt::from(n), BigInt::from(d))
}

/// The integer `n` as a rational.
///
/// # Examples
///
/// ```
/// use symplex::linprog::{q, qi};
/// assert_eq!(qi(3), q(6, 2));
/// ```
pub fn qi(n: i64) -> Q {
    Ratio::from_integer(BigInt::from(n))
}

fn invalid(operation: &'static str, reason: impl Into<String>) -> SymplexError {
    SymplexError::InvalidArgument {
        operation,
        reason: reason.into(),
    }
}

fn failed(operation: &'static str, reason: impl Into<String>) -> SymplexError {
    SymplexError::ComputationFailed {
        operation,
        reason: reason.into(),
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Problem description
// ═══════════════════════════════════════════════════════════════════════════

/// Optimisation direction.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum Objective {
    /// Minimise `cᵀx`.
    Minimize,
    /// Maximise `cᵀx`.
    Maximize,
}

/// Relation of a constraint row to its right-hand side.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Relation {
    Eq,
    Le,
    Ge,
}

#[derive(Clone, Debug)]
struct Constraint {
    row: Vec<Q>,
    rhs: Q,
    relation: Relation,
}

/// A linear program in builder form.  See the [module docs](self) for
/// the problem shape and the meaning of the results.
///
/// Every variable defaults to `0 ≤ xⱼ < ∞`; change that with
/// [`bounds`](Self::bounds) or [`free`](Self::free).  Constraint rows are
/// kept in insertion order, which is also the order of
/// [`LpSolution::duals`] and [`LpSolution::farkas`].
///
/// # Examples
///
/// ```
/// use symplex::linprog::{LpProblem, LpStatus, q, qi};
///
/// // min x + y  s.t.  x + 2y ≥ 1,  3x + y ≥ 1,  x, y ≥ 0   →  (1/5, 2/5)
/// let sol = LpProblem::minimize(vec![qi(1), qi(1)])
///     .ge(vec![qi(1), qi(2)], qi(1))
///     .ge(vec![qi(3), qi(1)], qi(1))
///     .solve()
///     .unwrap();
/// assert_eq!(sol.status, LpStatus::Optimal);
/// assert_eq!(sol.x, vec![q(1, 5), q(2, 5)]);
/// assert_eq!(sol.objective, Some(q(3, 5)));
/// ```
#[derive(Clone, Debug)]
pub struct LpProblem {
    objective: Objective,
    c: Vec<Q>,
    constraints: Vec<Constraint>,
    bounds: Vec<(Option<Q>, Option<Q>)>,
    /// First out-of-range variable index passed to `bounds`, reported by
    /// `solve`.
    bad_var: Option<usize>,
}

impl LpProblem {
    fn new(objective: Objective, c: Vec<Q>) -> Self {
        let n = c.len();
        LpProblem {
            objective,
            c,
            constraints: Vec::new(),
            bounds: vec![(Some(Q::zero()), None); n],
            bad_var: None,
        }
    }

    /// Start a minimisation of `cᵀx` over `x ∈ ℚ^{c.len()}`.
    pub fn minimize(c: Vec<Q>) -> Self {
        Self::new(Objective::Minimize, c)
    }

    /// Start a maximisation of `cᵀx` over `x ∈ ℚ^{c.len()}`.
    pub fn maximize(c: Vec<Q>) -> Self {
        Self::new(Objective::Maximize, c)
    }

    fn add(&mut self, row: Vec<Q>, rhs: Q, relation: Relation) {
        self.constraints.push(Constraint { row, rhs, relation });
    }

    /// Add the equality constraint `row·x = rhs`.
    pub fn eq(mut self, row: Vec<Q>, rhs: Q) -> Self {
        self.add(row, rhs, Relation::Eq);
        self
    }

    /// Add the inequality constraint `row·x ≤ rhs`.
    pub fn le(mut self, row: Vec<Q>, rhs: Q) -> Self {
        self.add(row, rhs, Relation::Le);
        self
    }

    /// Add the inequality constraint `row·x ≥ rhs`.
    pub fn ge(mut self, row: Vec<Q>, rhs: Q) -> Self {
        self.add(row, rhs, Relation::Ge);
        self
    }

    /// Set the bounds `lo ≤ x[var] ≤ hi` of one variable.  `None` means
    /// unbounded on that side; `(None, None)` makes the variable free.
    ///
    /// An out-of-range `var` is reported by [`solve`](Self::solve).
    pub fn bounds(mut self, var: usize, lo: Option<Q>, hi: Option<Q>) -> Self {
        match self.bounds.get_mut(var) {
            Some(slot) => *slot = (lo, hi),
            None => self.bad_var = self.bad_var.or(Some(var)),
        }
        self
    }

    /// Make `x[var]` free (`−∞ < x[var] < ∞`).
    pub fn free(self, var: usize) -> Self {
        self.bounds(var, None, None)
    }

    /// Number of decision variables (the length of `c`).
    pub fn num_vars(&self) -> usize {
        self.c.len()
    }

    /// Number of constraint rows added so far.
    pub fn num_constraints(&self) -> usize {
        self.constraints.len()
    }

    /// Solve the program exactly.
    ///
    /// The returned [`LpSolution`] reports [`LpStatus::Optimal`],
    /// [`LpStatus::Infeasible`] (with a Farkas certificate) or
    /// [`LpStatus::Unbounded`]; none of these is an error.
    ///
    /// # Errors
    ///
    /// * [`SymplexError::InvalidArgument`] for malformed input: no
    ///   variables, a constraint row whose length differs from `c.len()`,
    ///   or a bound set for a variable index `≥ c.len()`.
    /// * [`SymplexError::ComputationFailed`] if the pivot cap
    ///   (`10 000 + 50·(m + n)`) is exceeded — not expected in practice,
    ///   since Bland's rule rules out cycling.
    pub fn solve(&self) -> Result<LpSolution, SymplexError> {
        self.validate()?;
        solve_lp(self)
    }

    fn validate(&self) -> Result<(), SymplexError> {
        let n = self.c.len();
        if n == 0 {
            return Err(invalid(
                "linprog",
                "the objective must have at least one variable",
            ));
        }
        if let Some(bad) = self.bad_var {
            return Err(invalid(
                "linprog",
                format!("bounds were set for variable {bad} but there are only {n} variables"),
            ));
        }
        for (i, con) in self.constraints.iter().enumerate() {
            if con.row.len() != n {
                return Err(invalid(
                    "linprog",
                    format!(
                        "constraint {i} has {} coefficients but there are {n} variables",
                        con.row.len()
                    ),
                ));
            }
        }
        Ok(())
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Solution
// ═══════════════════════════════════════════════════════════════════════════

/// Outcome of a linear program.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum LpStatus {
    /// A finite optimum was found; `x`, `objective` and `duals` are set.
    Optimal,
    /// No point satisfies all constraints and bounds; `farkas` holds a
    /// certificate (see the [module docs](self)).
    Infeasible,
    /// The objective can be improved without limit over the feasible set.
    Unbounded,
}

/// Result of [`LpProblem::solve`] / [`linprog`].
///
/// See the [module docs](self) for the exact meaning and sign conventions
/// of `duals` and `farkas`.
#[derive(Clone, Debug)]
pub struct LpSolution {
    /// Optimal / infeasible / unbounded.
    pub status: LpStatus,
    /// The optimal point (empty unless `status == Optimal`).
    pub x: Vec<Q>,
    /// The optimal objective value `cᵀx` (`None` unless `Optimal`).
    pub objective: Option<Q>,
    /// Shadow prices, one per constraint in insertion order (empty unless
    /// `Optimal`): `yᵢ = ∂(optimal objective)/∂bᵢ`.
    pub duals: Vec<Q>,
    /// Farkas infeasibility certificate, one entry per constraint (`Some`
    /// only when `Infeasible` and the constraints — not the bounds alone —
    /// are contradictory).
    pub farkas: Option<Vec<Q>>,
}

impl LpSolution {
    /// `true` when `status == LpStatus::Optimal`.
    pub fn is_optimal(&self) -> bool {
        self.status == LpStatus::Optimal
    }

    /// The optimal point as exact rational expressions in `ctx`.
    pub fn x_ex(&self, ctx: &Context) -> Vec<Ex> {
        self.x.iter().map(|v| ctx.from_ratio(v.clone())).collect()
    }

    /// The shadow prices as exact rational expressions in `ctx`.
    pub fn duals_ex(&self, ctx: &Context) -> Vec<Ex> {
        self.duals
            .iter()
            .map(|v| ctx.from_ratio(v.clone()))
            .collect()
    }

    fn infeasible(farkas: Option<Vec<Q>>) -> Self {
        LpSolution {
            status: LpStatus::Infeasible,
            x: Vec::new(),
            objective: None,
            duals: Vec::new(),
            farkas,
        }
    }

    fn unbounded() -> Self {
        LpSolution {
            status: LpStatus::Unbounded,
            x: Vec::new(),
            objective: None,
            duals: Vec::new(),
            farkas: None,
        }
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Standard-form conversion
// ═══════════════════════════════════════════════════════════════════════════

/// How an original variable is represented in the standard form.
#[derive(Clone, Debug)]
enum VarMap {
    /// `x = lo + z[col]`, `z ≥ 0`.
    Shifted { col: usize, lo: Q },
    /// `x = hi − z[col]`, `z ≥ 0` (only an upper bound).
    Mirrored { col: usize, hi: Q },
    /// `x = z[pos] − z[neg]` (free variable).
    Split { pos: usize, neg: usize },
}

/// `min c·z  s.t.  A z = b, z ≥ 0, b ≥ 0`, plus the bookkeeping needed to
/// map results back to the caller's variables and constraints.
struct Standard {
    a: Vec<Vec<Q>>,
    b: Vec<Q>,
    c: Vec<Q>,
    /// Objective constant (min-form) introduced by variable shifts.
    constant: Q,
    var_map: Vec<VarMap>,
    /// `+1` / `−1`: whether original constraint `i` was negated to make
    /// its right-hand side non-negative.
    row_sign: Vec<Q>,
    /// Number of original constraints (leading rows of `a`).
    m_orig: usize,
}

/// Outcome of standardisation: either a standard-form problem or an
/// immediate verdict from the bounds.
enum Standardized {
    Ready(Standard),
    BoundsInfeasible,
}

fn standardize(p: &LpProblem) -> Standardized {
    let n = p.c.len();
    let m = p.constraints.len();
    let zero = Q::zero();
    let one = Q::one();

    // Variables.
    let mut var_map = Vec::with_capacity(n);
    let mut ncols = 0usize;
    // Rows for two-sided bounds `z ≤ hi − lo`, added after the constraints.
    let mut bound_rows: Vec<(usize, Q)> = Vec::new();
    for (lo, hi) in &p.bounds {
        match (lo, hi) {
            (Some(lo), Some(hi)) => {
                if hi < lo {
                    return Standardized::BoundsInfeasible;
                }
                var_map.push(VarMap::Shifted {
                    col: ncols,
                    lo: lo.clone(),
                });
                bound_rows.push((ncols, hi - lo));
                ncols += 1;
            }
            (Some(lo), None) => {
                var_map.push(VarMap::Shifted {
                    col: ncols,
                    lo: lo.clone(),
                });
                ncols += 1;
            }
            (None, Some(hi)) => {
                var_map.push(VarMap::Mirrored {
                    col: ncols,
                    hi: hi.clone(),
                });
                ncols += 1;
            }
            (None, None) => {
                var_map.push(VarMap::Split {
                    pos: ncols,
                    neg: ncols + 1,
                });
                ncols += 2;
            }
        }
    }
    let n_slack: usize = p
        .constraints
        .iter()
        .filter(|c| c.relation != Relation::Eq)
        .count();
    let total_cols = ncols + n_slack + bound_rows.len();
    let total_rows = m + bound_rows.len();

    // Objective in min form: (sign) · c, expressed in z.
    let sign = match p.objective {
        Objective::Minimize => one.clone(),
        Objective::Maximize => -one.clone(),
    };
    let mut c = vec![zero.clone(); total_cols];
    let mut constant = zero.clone();
    for (j, vm) in var_map.iter().enumerate() {
        let cj = &sign * &p.c[j];
        match vm {
            VarMap::Shifted { col, lo } => {
                c[*col] = cj.clone();
                constant += &cj * lo;
            }
            VarMap::Mirrored { col, hi } => {
                c[*col] = -cj.clone();
                constant += &cj * hi;
            }
            VarMap::Split { pos, neg } => {
                c[*pos] = cj.clone();
                c[*neg] = -cj;
            }
        }
    }

    // Constraint rows.
    let mut a = vec![vec![zero.clone(); total_cols]; total_rows];
    let mut b = vec![zero.clone(); total_rows];
    let mut row_sign = vec![one.clone(); m];
    let mut next_slack = ncols;
    for (i, con) in p.constraints.iter().enumerate() {
        let mut rhs = con.rhs.clone();
        for (j, vm) in var_map.iter().enumerate() {
            let aij = &con.row[j];
            if aij.is_zero() {
                continue;
            }
            match vm {
                VarMap::Shifted { col, lo } => {
                    a[i][*col] = aij.clone();
                    rhs -= aij * lo;
                }
                VarMap::Mirrored { col, hi } => {
                    a[i][*col] = -aij.clone();
                    rhs -= aij * hi;
                }
                VarMap::Split { pos, neg } => {
                    a[i][*pos] = aij.clone();
                    a[i][*neg] = -aij.clone();
                }
            }
        }
        match con.relation {
            Relation::Le => {
                a[i][next_slack] = one.clone();
                next_slack += 1;
            }
            Relation::Ge => {
                a[i][next_slack] = -one.clone();
                next_slack += 1;
            }
            Relation::Eq => {}
        }
        if rhs.is_negative() {
            row_sign[i] = -one.clone();
            for v in a[i].iter_mut() {
                if !v.is_zero() {
                    *v = -std::mem::take(v);
                }
            }
            rhs = -rhs;
        }
        b[i] = rhs;
    }
    for (k, (col, width)) in bound_rows.into_iter().enumerate() {
        let r = m + k;
        a[r][col] = one.clone();
        a[r][next_slack] = one.clone();
        next_slack += 1;
        b[r] = width;
    }
    debug_assert_eq!(next_slack, total_cols);

    Standardized::Ready(Standard {
        a,
        b,
        c,
        constant,
        var_map,
        row_sign,
        m_orig: m,
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// Simplex tableau
// ═══════════════════════════════════════════════════════════════════════════

/// Dense tableau for `min c·z s.t. A z = b, z ≥ 0` with one artificial
/// column per row.  Column layout: `[z columns | artificials | rhs]`.
///
/// **Integer pivoting.**  Every constraint row `i` is first multiplied by
/// `row_scale[i]`, the least common multiple of its denominators, so the
/// initial tableau is integral.  Pivots then follow Bareiss's fraction-free
/// Gauss–Jordan rule: all entries stay integers and share one common
/// denominator `d` (the current pivot, `±det B` of the scaled system).
/// The rational tableau entry is `rows[i][j] / d`; the reduced cost is
/// `obj[j] / (d · obj_scale)` where `obj_scale` clears the denominators of
/// the current cost vector.  Sign tests and the ratio test are integer
/// comparisons (cross-multiplied), so no gcd runs inside the pivot loop.
///
/// Scaling row `i` by `sᵢ` turns its artificial `aᵢ` into `a'ᵢ = sᵢ·aᵢ`.
/// Phase 1 therefore gives `a'ᵢ` the cost `1/sᵢ`, so that its objective is
/// still `Σ aᵢ`, every reduced cost is the same rational a `Ratio` tableau
/// of the unscaled system would hold, and the pivot sequence (Dantzig /
/// Bland choices, ratio tests, tie-breaks) is identical to it.
///
/// `obj` holds `d·obj_scale` times the reduced costs
/// `cⱼ − c_B B⁻¹ Aⱼ` and, in its last entry, `−(current objective value)`
/// at the same scale.  Because the artificial columns start as the
/// identity, at any time they hold `d·B⁻¹`, and the duals of the scaled
/// system are `y'ᵢ = c_{a'ᵢ} − obj[a'ᵢ]/(d·obj_scale)`; the duals of the
/// caller's rows are `yᵢ = sᵢ · y'ᵢ`.
///
/// **Hybrid arithmetic.**  The tableau is generic over its cell type
/// ([`Cell`]): it first runs on `i64` cells with `i128` intermediates, and
/// the first value that does not fit an `i64` aborts that attempt, after
/// which the whole problem is solved again on `BigInt` cells.  Every
/// decision (entering column, ratio test, tie-break) is a sign test or an
/// exact comparison of products, so the two runs take the same pivot path
/// and produce the same answer; the small-integer run merely never touches
/// the heap.  Certificate LPs (small polynomial coefficients) essentially
/// always stay in `i64`.
struct Tableau<T: Cell> {
    /// Row-major `m × width` integer tableau.
    rows: Vec<T>,
    width: usize,
    obj: Vec<T>,
    /// Common denominator of `rows` (sign may be negative after an
    /// artificial is driven out on a negative pivot).
    d: T,
    /// Positive integer clearing the denominators of the installed costs.
    obj_scale: T,
    /// Positive integer each constraint row was multiplied by.
    row_scale: Vec<T>,
    basis: Vec<usize>,
    m: usize,
    /// Number of genuine (non-artificial) columns.
    n: usize,
    pivots_done: usize,
    max_pivots: usize,
    /// Number of consecutive degenerate pivots (zero-length steps) so far.
    /// Dantzig's rule can only cycle through degenerate pivots, so after
    /// [`STALL_LIMIT`] of them in a row the entering rule switches to
    /// Bland's (which cannot cycle) until the next improving pivot resets
    /// the count.  Permanent Bland after the *first* degenerate pivot —
    /// the previous policy — made the highly degenerate certificate LPs
    /// (many zero right-hand sides) take tens of thousands of pivots.
    stall: usize,
}

/// Consecutive degenerate pivots tolerated under Dantzig's rule before
/// Bland's rule takes over for the rest of the stall.
const STALL_LIMIT: usize = 12;

/// A tableau cell: an integer type with the handful of exact operations
/// the fraction-free simplex needs.  Fallible operations return `None`
/// when the result does not fit; the `BigInt` implementation never fails.
trait Cell: Clone + PartialEq + Eq + Ord + fmt::Debug {
    fn cell_zero() -> Self;
    fn cell_one() -> Self;
    fn from_big(v: &BigInt) -> Option<Self>;
    fn to_big(&self) -> BigInt;
    fn is_zero(&self) -> bool;
    fn is_negative(&self) -> bool;
    /// `Ordering` of `self` against zero.
    fn signum(&self) -> std::cmp::Ordering;
    fn neg(&self) -> Option<Self>;
    fn mul(&self, o: &Self) -> Option<Self>;
    /// `(v · p − f · pr) / d`, the fraction-free pivot update; the division
    /// is exact.
    fn pivot_update(v: &Self, p: &Self, f: &Self, pr: &Self, d: &Self) -> Option<Self>;
    /// `(v · p) / d` (exact) — the update of a row with a zero pivot-column entry.
    fn rescale(v: &Self, p: &Self, d: &Self) -> Option<Self>;
    /// `self − f · r`, checked.
    fn sub_mul(&self, f: &Self, r: &Self) -> Option<Self>;
    /// `a · b` compared with `c · d`, exactly.
    fn cmp_products(a: &Self, b: &Self, c: &Self, d: &Self) -> std::cmp::Ordering;
}

impl Cell for BigInt {
    fn cell_zero() -> Self {
        <BigInt as Zero>::zero()
    }
    fn cell_one() -> Self {
        <BigInt as One>::one()
    }
    fn from_big(v: &BigInt) -> Option<Self> {
        Some(v.clone())
    }
    fn to_big(&self) -> BigInt {
        self.clone()
    }
    fn is_zero(&self) -> bool {
        Zero::is_zero(self)
    }
    fn is_negative(&self) -> bool {
        Signed::is_negative(self)
    }
    fn signum(&self) -> std::cmp::Ordering {
        match self.sign() {
            num_bigint::Sign::Minus => std::cmp::Ordering::Less,
            num_bigint::Sign::NoSign => std::cmp::Ordering::Equal,
            num_bigint::Sign::Plus => std::cmp::Ordering::Greater,
        }
    }
    fn neg(&self) -> Option<Self> {
        Some(-self)
    }
    fn mul(&self, o: &Self) -> Option<Self> {
        Some(self * o)
    }
    fn pivot_update(v: &Self, p: &Self, f: &Self, pr: &Self, d: &Self) -> Option<Self> {
        let t = if Zero::is_zero(pr) {
            v * p
        } else if Zero::is_zero(v) {
            -(f * pr)
        } else {
            v * p - f * pr
        };
        debug_assert!(
            Zero::is_zero(&(&t % d)),
            "integer pivoting: inexact division"
        );
        Some(t / d)
    }
    fn rescale(v: &Self, p: &Self, d: &Self) -> Option<Self> {
        let t = v * p;
        debug_assert!(
            Zero::is_zero(&(&t % d)),
            "integer pivoting: inexact division"
        );
        Some(t / d)
    }
    fn sub_mul(&self, f: &Self, r: &Self) -> Option<Self> {
        Some(self - f * r)
    }
    fn cmp_products(a: &Self, b: &Self, c: &Self, d: &Self) -> std::cmp::Ordering {
        (a * b).cmp(&(c * d))
    }
}

impl Cell for i64 {
    fn cell_zero() -> Self {
        0
    }
    fn cell_one() -> Self {
        1
    }
    fn from_big(v: &BigInt) -> Option<Self> {
        i64::try_from(v).ok()
    }
    fn to_big(&self) -> BigInt {
        BigInt::from(*self)
    }
    fn is_zero(&self) -> bool {
        *self == 0
    }
    fn is_negative(&self) -> bool {
        *self < 0
    }
    fn signum(&self) -> std::cmp::Ordering {
        self.cmp(&0)
    }
    fn neg(&self) -> Option<Self> {
        self.checked_neg()
    }
    fn mul(&self, o: &Self) -> Option<Self> {
        self.checked_mul(*o)
    }
    fn pivot_update(v: &Self, p: &Self, f: &Self, pr: &Self, d: &Self) -> Option<Self> {
        // Two i64 products and their difference always fit an i128.
        let t = i128::from(*v) * i128::from(*p) - i128::from(*f) * i128::from(*pr);
        debug_assert!(
            t % i128::from(*d) == 0,
            "integer pivoting: inexact division"
        );
        i64::try_from(t / i128::from(*d)).ok()
    }
    fn rescale(v: &Self, p: &Self, d: &Self) -> Option<Self> {
        let t = i128::from(*v) * i128::from(*p);
        debug_assert!(
            t % i128::from(*d) == 0,
            "integer pivoting: inexact division"
        );
        i64::try_from(t / i128::from(*d)).ok()
    }
    fn sub_mul(&self, f: &Self, r: &Self) -> Option<Self> {
        i64::try_from(i128::from(*self) - i128::from(*f) * i128::from(*r)).ok()
    }
    fn cmp_products(a: &Self, b: &Self, c: &Self, d: &Self) -> std::cmp::Ordering {
        (i128::from(*a) * i128::from(*b)).cmp(&(i128::from(*c) * i128::from(*d)))
    }
}

/// A signed 256-bit integer as `(hi, lo)` in two's complement, just enough
/// for the exact intermediates of `i128` cells: products of two `i128`s,
/// their differences, exact division by an `i128`, and comparison.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
struct I256 {
    hi: i128,
    lo: u128,
}

impl I256 {
    fn from_i128(v: i128) -> Self {
        I256 {
            hi: if v < 0 { -1 } else { 0 },
            lo: v as u128,
        }
    }

    /// Exact `a · b`.
    fn mul(a: i128, b: i128) -> Self {
        let neg = (a < 0) != (b < 0);
        let (ua, ub) = (a.unsigned_abs(), b.unsigned_abs());
        // 128×128 → 256 by 64-bit limbs.
        let (a0, a1) = (ua as u64 as u128, ua >> 64);
        let (b0, b1) = (ub as u64 as u128, ub >> 64);
        let p00 = a0 * b0;
        let p01 = a0 * b1;
        let p10 = a1 * b0;
        let p11 = a1 * b1;
        let mid = (p00 >> 64) + (p01 as u64 as u128) + (p10 as u64 as u128);
        let lo = (p00 as u64 as u128) | (mid << 64);
        let hi = p11 + (p01 >> 64) + (p10 >> 64) + (mid >> 64);
        let mag = I256 { hi: hi as i128, lo };
        if neg { mag.neg() } else { mag }
    }

    fn neg(self) -> Self {
        let lo = (!self.lo).wrapping_add(1);
        let hi = (!self.hi).wrapping_add(u128::from(lo == 0) as i128);
        I256 { hi, lo }
    }

    fn sub(self, o: Self) -> Self {
        let (lo, borrow) = self.lo.overflowing_sub(o.lo);
        let hi = self.hi.wrapping_sub(o.hi).wrapping_sub(i128::from(borrow));
        I256 { hi, lo }
    }

    fn is_negative(self) -> bool {
        self.hi < 0
    }

    /// The value as an `i128`, if it fits.
    fn to_i128(self) -> Option<i128> {
        let lo = self.lo as i128;
        let fits = (self.hi == 0 && lo >= 0) || (self.hi == -1 && lo < 0);
        fits.then_some(lo)
    }

    /// Exact division of a non-negative `self` by a positive `d` (the
    /// remainder is known to be zero); `None` if the quotient does not fit
    /// a `u128`.
    ///
    /// Jebelean's exact division: write `d = 2^k · d'` with `d'` odd, shift
    /// the dividend right by `k` (exact), and multiply its low 128 bits by
    /// the inverse of `d'` modulo `2^128` — which is the quotient whenever
    /// the quotient fits, i.e. whenever the shifted high half is below
    /// `d'`.  A handful of wrapping multiplications instead of a long
    /// division.
    fn div_exact_unsigned(self, d: u128) -> Option<u128> {
        debug_assert!(!self.is_negative());
        let hi = self.hi as u128;
        if hi == 0 {
            return Some(self.lo / d);
        }
        let k = d.trailing_zeros();
        let d_odd = d >> k;
        let lo = if k == 0 {
            self.lo
        } else {
            (self.lo >> k) | (hi << (128 - k))
        };
        let hi = hi >> k;
        if hi >= d_odd {
            return None;
        }
        Some(lo.wrapping_mul(inverse_mod_2_128(d_odd)))
    }

    /// Exact `self / d` for any signs, `None` if the quotient does not fit
    /// an `i128`.
    fn div_exact(self, d: i128) -> Option<i128> {
        let neg = self.is_negative() != (d < 0);
        let mag = if self.is_negative() { self.neg() } else { self };
        let q = mag.div_exact_unsigned(d.unsigned_abs())?;
        if neg {
            // −2^127 is representable; q up to 2^127 is allowed then.
            if q > (1u128 << 127) {
                return None;
            }
            Some((q as i128).wrapping_neg())
        } else {
            i128::try_from(q).ok()
        }
    }
}

/// The inverse of an odd `d` modulo `2^128`, by Newton's iteration
/// `x ← x·(2 − d·x)`, which doubles the number of correct low bits each
/// step (an odd `d` is its own inverse modulo 8, so three bits to start).
fn inverse_mod_2_128(d: u128) -> u128 {
    debug_assert!(d & 1 == 1);
    let mut x = d;
    for _ in 0..6 {
        x = x.wrapping_mul(2u128.wrapping_sub(d.wrapping_mul(x)));
    }
    debug_assert_eq!(d.wrapping_mul(x), 1);
    x
}

impl PartialOrd for I256 {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for I256 {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        self.hi.cmp(&other.hi).then(self.lo.cmp(&other.lo))
    }
}

impl Cell for i128 {
    fn cell_zero() -> Self {
        0
    }
    fn cell_one() -> Self {
        1
    }
    fn from_big(v: &BigInt) -> Option<Self> {
        i128::try_from(v).ok()
    }
    fn to_big(&self) -> BigInt {
        BigInt::from(*self)
    }
    fn is_zero(&self) -> bool {
        *self == 0
    }
    fn is_negative(&self) -> bool {
        *self < 0
    }
    fn signum(&self) -> std::cmp::Ordering {
        self.cmp(&0)
    }
    fn neg(&self) -> Option<Self> {
        self.checked_neg()
    }
    fn mul(&self, o: &Self) -> Option<Self> {
        self.checked_mul(*o)
    }
    fn pivot_update(v: &Self, p: &Self, f: &Self, pr: &Self, d: &Self) -> Option<Self> {
        let t = I256::mul(*v, *p).sub(I256::mul(*f, *pr));
        t.div_exact(*d)
    }
    fn rescale(v: &Self, p: &Self, d: &Self) -> Option<Self> {
        I256::mul(*v, *p).div_exact(*d)
    }
    fn sub_mul(&self, f: &Self, r: &Self) -> Option<Self> {
        I256::from_i128(*self).sub(I256::mul(*f, *r)).to_i128()
    }
    fn cmp_products(a: &Self, b: &Self, c: &Self, d: &Self) -> std::cmp::Ordering {
        I256::mul(*a, *b).cmp(&I256::mul(*c, *d))
    }
}

enum Step {
    Optimal,
    Unbounded,
}

/// Why a tableau run stopped early: a value did not fit the cell type
/// (retry with a wider one), or the pivot cap was hit.
enum Halt {
    Overflow,
    Error(SymplexError),
}

impl From<SymplexError> for Halt {
    fn from(e: SymplexError) -> Self {
        Halt::Error(e)
    }
}

/// Least common multiple of the denominators of `values`.
fn denominator_lcm<'a>(values: impl Iterator<Item = &'a Q>) -> BigInt {
    values.fold(<BigInt as One>::one(), |l, q| l.lcm(q.denom()))
}

impl<T: Cell> Tableau<T> {
    /// Build the initial tableau; `None` if an entry does not fit `T`.
    fn new(std: &Standard) -> Option<Self> {
        let m = std.a.len();
        let n = std.c.len();
        let width = n + m + 1;
        let mut rows = Vec::with_capacity(m * width);
        let mut row_scale = Vec::with_capacity(m);
        for i in 0..m {
            let s = denominator_lcm(std.a[i].iter().chain(std::iter::once(&std.b[i])));
            let scaled = |q: &Q| T::from_big(&(q.numer() * (&s / q.denom())));
            for q in &std.a[i] {
                rows.push(scaled(q)?);
            }
            for k in 0..m {
                rows.push(if k == i {
                    T::cell_one()
                } else {
                    T::cell_zero()
                });
            }
            rows.push(scaled(&std.b[i])?);
            row_scale.push(T::from_big(&s)?);
        }
        let basis: Vec<usize> = (0..m).map(|i| n + i).collect();
        Some(Tableau {
            rows,
            width,
            obj: vec![T::cell_zero(); width],
            d: T::cell_one(),
            obj_scale: T::cell_one(),
            row_scale,
            basis,
            m,
            n,
            pivots_done: 0,
            max_pivots: 10_000 + 50 * (m + n),
            stall: 0,
        })
    }

    #[inline]
    fn rhs_col(&self) -> usize {
        self.n + self.m
    }

    #[inline]
    fn row(&self, i: usize) -> &[T] {
        &self.rows[i * self.width..(i + 1) * self.width]
    }

    /// Sign of the rational value `z / d`.
    #[inline]
    fn sign_of(&self, z: &T) -> std::cmp::Ordering {
        use std::cmp::Ordering::*;
        match (z.signum(), self.d.signum()) {
            (Equal, _) => Equal,
            (a, b) if a == b => Greater,
            _ => Less,
        }
    }

    /// Rational value of a tableau entry.
    #[inline]
    fn value(&self, z: &T) -> Q {
        Ratio::new(z.to_big(), self.d.to_big())
    }

    /// Install a new objective (`costs` over all `n + m` columns) and price
    /// out the current basis so that `obj` holds true reduced costs (times
    /// `d · obj_scale`).  `None` on overflow.
    fn set_objective(&mut self, costs: &[Q]) -> Option<()> {
        let width = self.width;
        let s = denominator_lcm(costs.iter());
        let int_cost = |q: &Q| T::from_big(&(q.numer() * (&s / q.denom())));
        let mut obj = vec![T::cell_zero(); width];
        for (o, c) in obj.iter_mut().zip(costs) {
            if !Zero::is_zero(c) {
                *o = int_cost(c)?.mul(&self.d)?;
            }
        }
        for (i, &k) in self.basis.iter().enumerate() {
            let f = int_cost(&costs[k])?;
            if f.is_zero() {
                continue;
            }
            let row = self.row(i);
            for (o, r) in obj.iter_mut().zip(row) {
                if !r.is_zero() {
                    *o = o.sub_mul(&f, r)?;
                }
            }
        }
        self.obj = obj;
        self.obj_scale = T::from_big(&s)?;
        Some(())
    }

    /// Fraction-free pivot on `(r, s)`: every other row `i` (and the
    /// objective row) becomes `(row_i·p − row_i[s]·row_r) / d`, the pivot
    /// row is unchanged, and `d ← p`.  `None` on overflow (the tableau is
    /// then unusable; the caller restarts with wider cells).
    fn pivot(&mut self, r: usize, s: usize) -> Option<()> {
        let w = self.width;
        let prow: Vec<T> = self.rows[r * w..(r + 1) * w]
            .iter_mut()
            .map(|v| std::mem::replace(v, T::cell_zero()))
            .collect();
        let p = prow[s].clone();
        let d = std::mem::replace(&mut self.d, p.clone());
        let update_row = |row: &mut [T]| -> Option<()> {
            let f = std::mem::replace(&mut row[s], T::cell_zero());
            if f.is_zero() {
                for v in row.iter_mut() {
                    if !v.is_zero() {
                        *v = T::rescale(v, &p, &d)?;
                    }
                }
                return Some(());
            }
            for (j, (v, pr)) in row.iter_mut().zip(prow.iter()).enumerate() {
                if j == s {
                    continue;
                }
                if v.is_zero() && pr.is_zero() {
                    continue;
                }
                *v = T::pivot_update(v, &p, &f, pr, &d)?;
            }
            Some(())
        };
        for i in 0..self.m {
            if i == r {
                continue;
            }
            update_row(&mut self.rows[i * w..(i + 1) * w])?;
        }
        update_row(&mut self.obj)?;
        self.rows[r * w..(r + 1) * w]
            .iter_mut()
            .zip(prow)
            .for_each(|(slot, v)| *slot = v);
        self.basis[r] = s;
        self.pivots_done += 1;
        Some(())
    }

    /// Run the simplex method on the current objective, considering only
    /// columns `< limit` as candidates to enter the basis.
    fn run(&mut self, limit: usize) -> Result<Step, Halt> {
        use std::cmp::Ordering;
        loop {
            if self.pivots_done >= self.max_pivots {
                return Err(failed(
                    "linprog",
                    format!(
                        "pivot cap of {} exceeded; the problem is degenerate beyond \
                         what the solver handles",
                        self.max_pivots
                    ),
                )
                .into());
            }
            // Entering column: Dantzig (most negative reduced cost, lowest
            // index on ties) or Bland (lowest index with negative cost).
            // Reduced costs share the positive factor |d|·obj_scale, so the
            // integer entries compare like the rationals once oriented by
            // the sign of d.  The scaled artificial a'ᵢ = sᵢ·aᵢ has reduced
            // cost rc(aᵢ)/sᵢ; multiplying by sᵢ compares the reduced cost of
            // the *unscaled* artificial, so Dantzig's choice is exactly the
            // one a `Ratio` tableau of the caller's system would make.
            let neg = self.d.is_negative();
            let oriented = |z: &T, j: usize| -> Option<T> {
                let v = if neg { z.neg()? } else { z.clone() };
                if j >= self.n {
                    v.mul(&self.row_scale[j - self.n])
                } else {
                    Some(v)
                }
            };
            let mut entering: Option<(usize, T)> = None;
            for j in 0..limit {
                if self.sign_of(&self.obj[j]) != Ordering::Less {
                    continue;
                }
                let v = oriented(&self.obj[j], j).ok_or(Halt::Overflow)?;
                match &entering {
                    Some((_, best)) if v >= *best => {}
                    _ => entering = Some((j, v)),
                }
                if self.stall >= STALL_LIMIT {
                    break;
                }
            }
            let Some((s, _)) = entering else {
                return Ok(Step::Optimal);
            };
            // Leaving row: minimum ratio rhsᵢ/aᵢₛ (d cancels), ties broken by
            // smallest basic index.  Every eligible aᵢₛ has the sign of d, so
            // rhsᵢ/aᵢₛ < rhsₗ/aₗₛ  ⇔  rhsᵢ·aₗₛ < rhsₗ·aᵢₛ  (cross-multiplication
            // by a positive product; no gcd needed).
            let rhs = self.rhs_col();
            let mut leaving: Option<usize> = None;
            for i in 0..self.m {
                let a = &self.row(i)[s];
                if self.sign_of(a) != Ordering::Greater {
                    continue;
                }
                let better = match leaving {
                    None => true,
                    Some(l) => {
                        let (ri, rl) = (&self.row(i)[rhs], &self.row(l)[rhs]);
                        let al = &self.row(l)[s];
                        match T::cmp_products(ri, al, rl, a) {
                            Ordering::Less => true,
                            Ordering::Equal => self.basis[i] < self.basis[l],
                            Ordering::Greater => false,
                        }
                    }
                };
                if better {
                    leaving = Some(i);
                }
            }
            let Some(r) = leaving else {
                return Ok(Step::Unbounded);
            };
            if self.row(r)[rhs].is_zero() {
                self.stall += 1;
            } else {
                self.stall = 0;
            }
            self.pivot(r, s).ok_or(Halt::Overflow)?;
        }
    }

    /// Value of the current basic solution over the genuine columns.
    fn solution(&self) -> Vec<Q> {
        let rhs = self.rhs_col();
        let mut z = vec![Q::zero(); self.n];
        for (i, &k) in self.basis.iter().enumerate() {
            if k < self.n {
                z[k] = self.value(&self.row(i)[rhs]);
            }
        }
        z
    }

    /// `−(current objective value)` — the last entry of the objective row.
    fn neg_objective(&self) -> Q {
        Ratio::new(
            self.obj[self.rhs_col()].to_big(),
            self.d.to_big() * self.obj_scale.to_big(),
        )
    }

    /// Pivot artificial variables out of the basis wherever possible.
    /// Rows whose genuine part is entirely zero are redundant; their
    /// artificial stays basic at level zero and never moves again.
    fn drive_out_artificials(&mut self) -> Option<()> {
        for r in 0..self.m {
            if self.basis[r] < self.n {
                continue;
            }
            if let Some(s) = (0..self.n).find(|&j| !self.row(r)[j].is_zero()) {
                self.pivot(r, s)?;
            }
        }
        Some(())
    }

    /// Phase-1 cost vector over all `n + m` columns: `0` for genuine
    /// columns and `1/sᵢ` for the (scaled) artificial of row `i`, so that
    /// the phase-1 objective is `Σ aᵢ` of the unscaled system.
    fn phase1_costs(&self) -> Vec<Q> {
        let mut costs = vec![Q::zero(); self.n + self.m];
        for (c, s) in costs[self.n..].iter_mut().zip(&self.row_scale) {
            *c = Ratio::new(<BigInt as One>::one(), s.to_big());
        }
        costs
    }

    /// Current duals of the caller's (unscaled) rows: `yᵢ = sᵢ · (c_{a'ᵢ} −
    /// reduced cost of a'ᵢ)`, where the artificial's cost is `1/sᵢ` in
    /// phase 1 and `0` in phase 2.
    fn duals(&self, phase1: bool) -> Vec<Q> {
        let denom = self.d.to_big() * self.obj_scale.to_big();
        (0..self.m)
            .map(|i| {
                let s = self.row_scale[i].to_big();
                let reduced = Ratio::new(self.obj[self.n + i].to_big(), denom.clone());
                let scaled_reduced = reduced * Ratio::from_integer(s);
                if phase1 {
                    Q::one() - scaled_reduced
                } else {
                    -scaled_reduced
                }
            })
            .collect()
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Driver
// ═══════════════════════════════════════════════════════════════════════════

fn solve_lp(p: &LpProblem) -> Result<LpSolution, SymplexError> {
    let sf = match standardize(p) {
        Standardized::Ready(s) => s,
        Standardized::BoundsInfeasible => return Ok(LpSolution::infeasible(None)),
    };
    // Fixed-width cells first (i64, then i128 with exact 256-bit
    // intermediates); the same algorithm on BigInt cells if a value
    // outgrows them.  Fraction-free entries are minors of the scaled
    // system, so a 20-row tableau typically peaks around 70–100 bits.
    match solve_standard::<i64>(p, &sf) {
        Ok(sol) => return Ok(sol),
        Err(Halt::Error(e)) => return Err(e),
        Err(Halt::Overflow) => {}
    }
    match solve_standard::<i128>(p, &sf) {
        Ok(sol) => return Ok(sol),
        Err(Halt::Error(e)) => return Err(e),
        Err(Halt::Overflow) => {}
    }
    tracing::debug!(target: "symplex::linprog", rows = sf.a.len(), cols = sf.c.len(), "i128 tableau overflowed; solving on BigInt");
    match solve_standard::<BigInt>(p, &sf) {
        Ok(sol) => Ok(sol),
        Err(Halt::Error(e)) => Err(e),
        Err(Halt::Overflow) => Err(failed("linprog", "internal: BigInt tableau overflowed")),
    }
}

fn solve_standard<T: Cell>(p: &LpProblem, sf: &Standard) -> Result<LpSolution, Halt> {
    let mut t = Tableau::<T>::new(sf).ok_or(Halt::Overflow)?;
    let m = t.m;
    let n = t.n;

    // Phase 1: minimise the sum of the artificials.
    let phase1 = t.phase1_costs();
    t.set_objective(&phase1).ok_or(Halt::Overflow)?;
    match t.run(n + m)? {
        Step::Optimal => {}
        Step::Unbounded => {
            // Cannot happen: the phase-1 objective is bounded below by 0.
            return Err(failed("linprog", "phase 1 reported unbounded").into());
        }
    }
    let infeasibility = -t.neg_objective();
    if infeasibility.is_positive() {
        // Farkas certificate from the phase-1 duals, mapped back to the
        // caller's rows (bound rows are absorbed into the box infimum).
        let y_std = t.duals(true);
        let farkas: Vec<Q> = (0..sf.m_orig)
            .map(|i| -(&sf.row_sign[i] * &y_std[i]))
            .collect();
        return Ok(LpSolution::infeasible(Some(farkas)));
    }

    // Phase 2.
    t.drive_out_artificials().ok_or(Halt::Overflow)?;
    let mut phase2 = vec![Q::zero(); n + m];
    phase2[..n].clone_from_slice(&sf.c);
    t.set_objective(&phase2).ok_or(Halt::Overflow)?;
    match t.run(n)? {
        Step::Optimal => {}
        Step::Unbounded => return Ok(LpSolution::unbounded()),
    }

    if tracing::enabled!(target: "symplex::linprog::growth", tracing::Level::TRACE) {
        let bits = t
            .rows
            .iter()
            .chain(t.obj.iter())
            .chain(std::iter::once(&t.d))
            .map(|c| c.to_big().bits())
            .max()
            .unwrap_or(0);
        tracing::trace!(target: "symplex::linprog::growth", rows = t.m, cols = t.n, pivots = t.pivots_done, max_bits = bits, "final tableau");
    }
    let z = t.solution();
    let x: Vec<Q> = sf
        .var_map
        .iter()
        .map(|vm| match vm {
            VarMap::Shifted { col, lo } => lo + &z[*col],
            VarMap::Mirrored { col, hi } => hi - &z[*col],
            VarMap::Split { pos, neg } => &z[*pos] - &z[*neg],
        })
        .collect();
    let objective: Q = p.c.iter().zip(x.iter()).map(|(c, v)| c * v).sum();
    // The tableau's objective (min form, without the shift constant) must
    // agree with the value recomputed from x.
    debug_assert_eq!(
        {
            let min_form: Q =
                sf.c.iter().zip(z.iter()).map(|(c, v)| c * v).sum::<Q>() + &sf.constant;
            match p.objective {
                Objective::Minimize => min_form,
                Objective::Maximize => -min_form,
            }
        },
        objective
    );

    let y_std = t.duals(false);
    let dir = match p.objective {
        Objective::Minimize => Q::one(),
        Objective::Maximize => -Q::one(),
    };
    let duals: Vec<Q> = (0..sf.m_orig)
        .map(|i| &dir * &(&sf.row_sign[i] * &y_std[i]))
        .collect();

    Ok(LpSolution {
        status: LpStatus::Optimal,
        x,
        objective: Some(objective),
        duals,
        farkas: None,
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// Convenience wrappers
// ═══════════════════════════════════════════════════════════════════════════

/// SciPy-shaped entry point: minimise `cᵀx` subject to `A_ub x ≤ b_ub`,
/// `A_eq x = b_eq` and `bounds[j].0 ≤ xⱼ ≤ bounds[j].1`.
///
/// `bounds` may be empty (every variable then defaults to `0 ≤ xⱼ < ∞`) or
/// have one `(lo, hi)` pair per variable, `None` meaning unbounded on that
/// side.  Constraints are numbered with the `≤` rows first, then the `=`
/// rows — that is the order of [`LpSolution::duals`] / [`farkas`](LpSolution::farkas).
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for malformed input (empty `c`, a row
/// of the wrong length, `b` of the wrong length, `bounds` of the wrong
/// length).
///
/// # Examples
///
/// ```
/// use symplex::linprog::{linprog, LpStatus, q, qi};
///
/// // min −x − y   s.t.  x + 2y ≤ 4,  3x + y ≤ 6,  x, y ≥ 0
/// let sol = linprog(
///     &[qi(-1), qi(-1)],
///     &[vec![qi(1), qi(2)], vec![qi(3), qi(1)]],
///     &[qi(4), qi(6)],
///     &[],
///     &[],
///     &[],
/// )
/// .unwrap();
/// assert_eq!(sol.status, LpStatus::Optimal);
/// assert_eq!(sol.x, vec![q(8, 5), q(6, 5)]);
/// assert_eq!(sol.objective, Some(q(-14, 5)));
/// ```
pub fn linprog(
    c: &[Q],
    a_ub: &[Vec<Q>],
    b_ub: &[Q],
    a_eq: &[Vec<Q>],
    b_eq: &[Q],
    bounds: &[(Option<Q>, Option<Q>)],
) -> Result<LpSolution, SymplexError> {
    if a_ub.len() != b_ub.len() {
        return Err(invalid(
            "linprog",
            format!(
                "A_ub has {} rows but b_ub has {} entries",
                a_ub.len(),
                b_ub.len()
            ),
        ));
    }
    if a_eq.len() != b_eq.len() {
        return Err(invalid(
            "linprog",
            format!(
                "A_eq has {} rows but b_eq has {} entries",
                a_eq.len(),
                b_eq.len()
            ),
        ));
    }
    if !bounds.is_empty() && bounds.len() != c.len() {
        return Err(invalid(
            "linprog",
            format!(
                "bounds has {} entries but there are {} variables",
                bounds.len(),
                c.len()
            ),
        ));
    }
    let mut p = LpProblem::minimize(c.to_vec());
    for (row, rhs) in a_ub.iter().zip(b_ub) {
        p = p.le(row.clone(), rhs.clone());
    }
    for (row, rhs) in a_eq.iter().zip(b_eq) {
        p = p.eq(row.clone(), rhs.clone());
    }
    for (j, (lo, hi)) in bounds.iter().enumerate() {
        p = p.bounds(j, lo.clone(), hi.clone());
    }
    p.solve()
}

/// One-line human-readable summary: status, point, objective and duals.
///
/// ```
/// use symplex::linprog::{LpProblem, qi};
///
/// let sol = LpProblem::maximize(vec![qi(3), qi(2)])
///     .le(vec![qi(1), qi(1)], qi(4))
///     .le(vec![qi(1), qi(3)], qi(6))
///     .solve()
///     .unwrap();
/// assert_eq!(sol.to_string(), "Optimal: x = (4, 0), objective = 12, duals = (3, 0)");
/// ```
impl fmt::Display for LpSolution {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let list = |v: &[Q]| {
            let parts: Vec<String> = v.iter().map(ToString::to_string).collect();
            format!("({})", parts.join(", "))
        };
        match self.status {
            LpStatus::Optimal => {
                write!(f, "Optimal: x = {}", list(&self.x))?;
                if let Some(obj) = &self.objective {
                    write!(f, ", objective = {obj}")?;
                }
                write!(f, ", duals = {}", list(&self.duals))
            }
            LpStatus::Infeasible => match &self.farkas {
                Some(y) => write!(f, "Infeasible: Farkas certificate y = {}", list(y)),
                None => write!(f, "Infeasible: contradictory bounds"),
            },
            LpStatus::Unbounded => write!(f, "Unbounded"),
        }
    }
}

/// Outcome of an exact feasibility question (see [`feasible_nonneg_certified`]
/// and [`nonneg_combination`]).
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum Feasibility {
    /// A witness `x ≥ 0` with `A x = b`.
    Feasible(Vec<Q>),
    /// No such `x`; `farkas` is a vector `y` with `Aᵀy ≥ 0` and `yᵀb < 0`
    /// (one entry per equation), which proves it.
    Infeasible {
        /// The Farkas certificate, `None` only if the solver could not
        /// produce one (never the case for pure equality systems).
        farkas: Option<Vec<Q>>,
    },
}

impl Feasibility {
    /// The witness, if feasible.
    pub fn witness(&self) -> Option<&[Q]> {
        match self {
            Feasibility::Feasible(x) => Some(x),
            Feasibility::Infeasible { .. } => None,
        }
    }

    /// `true` for [`Feasibility::Feasible`].
    pub fn is_feasible(&self) -> bool {
        matches!(self, Feasibility::Feasible(_))
    }
}

/// Exact feasibility of `A x = b, x ≥ 0`: `Some(x)` with a solution, or
/// `None` if none exists.  Use [`feasible_nonneg_certified`] when you also
/// want the Farkas certificate for the infeasible case.
///
/// This is the question behind many certificate searches (Farkas,
/// Positivstellensatz-style combinations, Carathéodory decompositions);
/// floating-point solvers answer it only up to a tolerance.
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if `a_eq` is empty, jagged, or its
/// length differs from `b_eq`.
///
/// # Examples
///
/// ```
/// use symplex::linprog::{feasible_nonneg, q, qi};
///
/// // x/3 + y/7 = 1,  x − y = 0  →  x = y = 21/10
/// let x = feasible_nonneg(
///     &[vec![q(1, 3), q(1, 7)], vec![qi(1), qi(-1)]],
///     &[qi(1), qi(0)],
/// )
/// .unwrap()
/// .unwrap();
/// assert_eq!(x, vec![q(21, 10), q(21, 10)]);
///
/// // x + y = −1 has no non-negative solution.
/// assert!(feasible_nonneg(&[vec![qi(1), qi(1)]], &[qi(-1)]).unwrap().is_none());
/// ```
pub fn feasible_nonneg(a_eq: &[Vec<Q>], b_eq: &[Q]) -> Result<Option<Vec<Q>>, SymplexError> {
    let Some(first) = a_eq.first() else {
        return Err(invalid(
            "feasible_nonneg",
            "need at least one equation to determine the number of variables",
        ));
    };
    let n = first.len();
    if n == 0 {
        return Err(invalid(
            "feasible_nonneg",
            "equations must have at least one variable",
        ));
    }
    let sol = linprog(&vec![Q::zero(); n], &[], &[], a_eq, b_eq, &[])?;
    Ok(match sol.status {
        LpStatus::Optimal => Some(sol.x),
        LpStatus::Infeasible => None,
        // A zero objective cannot be unbounded.
        LpStatus::Unbounded => None,
    })
}

/// Like [`feasible_nonneg`], but an infeasible system comes back with its
/// Farkas certificate: a `y` (one entry per equation) with `Aᵀy ≥ 0`
/// component-wise and `yᵀb < 0`, so `x ≥ 0 ⇒ yᵀ(Ax) ≥ 0 > yᵀb`.
///
/// # Errors
///
/// As [`feasible_nonneg`].
///
/// # Examples
///
/// ```
/// use symplex::linprog::{feasible_nonneg_certified, Feasibility, qi};
///
/// // x + y = 1  and  x + y = 2 cannot both hold.
/// let a = [vec![qi(1), qi(1)], vec![qi(1), qi(1)]];
/// match feasible_nonneg_certified(&a, &[qi(1), qi(2)]).unwrap() {
///     Feasibility::Infeasible { farkas: Some(y) } => {
///         // Aᵀy = (y₀ + y₁, y₀ + y₁) ≥ 0  and  y₀ + 2·y₁ < 0
///         let g = &y[0] + &y[1];
///         assert!(g >= qi(0));
///         assert!(&y[0] + &y[1] * qi(2) < qi(0));
///     }
///     other => panic!("expected a certificate, got {other:?}"),
/// }
/// ```
pub fn feasible_nonneg_certified(a_eq: &[Vec<Q>], b_eq: &[Q]) -> Result<Feasibility, SymplexError> {
    let Some(first) = a_eq.first() else {
        return Err(invalid(
            "feasible_nonneg_certified",
            "need at least one equation to determine the number of variables",
        ));
    };
    let n = first.len();
    if n == 0 {
        return Err(invalid(
            "feasible_nonneg_certified",
            "equations must have at least one variable",
        ));
    }
    let sol = linprog(&vec![Q::zero(); n], &[], &[], a_eq, b_eq, &[])?;
    Ok(match sol.status {
        LpStatus::Optimal => Feasibility::Feasible(sol.x),
        LpStatus::Infeasible => Feasibility::Infeasible { farkas: sol.farkas },
        // A zero objective cannot be unbounded.
        LpStatus::Unbounded => Feasibility::Infeasible { farkas: None },
    })
}

/// Is `target` a non-negative combination `Σ λⱼ vⱼ` of the `vectors`?
///
/// The cone-membership form of [`feasible_nonneg_certified`]: the vectors
/// are the *columns* of `A` (each `vectors[j]` has one entry per coordinate
/// of `target`), and the answer is the coefficient vector `λ ≥ 0` or a
/// Farkas certificate `y` with `y·vⱼ ≥ 0` for every `j` and `y·target < 0`
/// (a hyperplane separating `target` from the cone).
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if `vectors` is empty or some vector
/// has a different length than `target`.
///
/// # Examples
///
/// ```
/// use symplex::linprog::{nonneg_combination, Feasibility, q, qi};
///
/// // (1, 1) = ½·(2, 0) + 1·(0, 1)
/// let cone = [vec![qi(2), qi(0)], vec![qi(0), qi(1)]];
/// assert_eq!(
///     nonneg_combination(&cone, &[qi(1), qi(1)]).unwrap(),
///     Feasibility::Feasible(vec![q(1, 2), qi(1)])
/// );
/// // (−1, 1) is outside the cone spanned by (2, 0) and (0, 1).
/// assert!(!nonneg_combination(&cone, &[qi(-1), qi(1)]).unwrap().is_feasible());
/// ```
pub fn nonneg_combination(vectors: &[Vec<Q>], target: &[Q]) -> Result<Feasibility, SymplexError> {
    if vectors.is_empty() {
        return Err(invalid("nonneg_combination", "need at least one vector"));
    }
    let m = target.len();
    if m == 0 {
        return Err(invalid(
            "nonneg_combination",
            "target must have at least one coordinate",
        ));
    }
    if let Some((j, v)) = vectors.iter().enumerate().find(|(_, v)| v.len() != m) {
        return Err(invalid(
            "nonneg_combination",
            format!(
                "vector {j} has {} coordinates but the target has {m}",
                v.len()
            ),
        ));
    }
    // Rows of A are coordinates, columns are the vectors.
    let a_eq: Vec<Vec<Q>> = (0..m)
        .map(|i| vectors.iter().map(|v| v[i].clone()).collect())
        .collect();
    feasible_nonneg_certified(&a_eq, target).map_err(|e| match e {
        SymplexError::InvalidArgument { reason, .. } => invalid("nonneg_combination", reason),
        other => other,
    })
}

/// Rational value of every entry of `m` (after constant folding), or an
/// `InvalidArgument` naming the offending entry.
fn matrix_to_q(m: &Matrix, what: &str) -> Result<Vec<Vec<Q>>, SymplexError> {
    if let Some(rows) = m.to_rational_rows() {
        return Ok(rows);
    }
    let evaluated = m.eval();
    evaluated.to_rational_rows().ok_or_else(|| {
        let bad = evaluated
            .iter()
            .find(|e| e.as_rational().is_none())
            .map(|e| e.to_string())
            .unwrap_or_default();
        invalid(
            "linprog_matrix",
            format!("{what} must contain only numeric literals; found `{bad}`"),
        )
    })
}

/// A column or row vector as a flat list.
fn matrix_to_vec(m: &Matrix, what: &str) -> Result<Vec<Q>, SymplexError> {
    if m.ncols() != 1 && m.nrows() != 1 {
        return Err(invalid(
            "linprog_matrix",
            format!(
                "{what} must be a row or column vector, got {}×{}",
                m.nrows(),
                m.ncols()
            ),
        ));
    }
    // Either one row or one entry per row: flattening gives the entries in
    // order in both cases.
    Ok(matrix_to_q(m, what)?.into_iter().flatten().collect())
}

/// Append the rows of `A x (relation) b` to `p`, checking shapes.
fn add_matrix_block(
    p: &mut LpProblem,
    a: Option<&Matrix>,
    b: Option<&Matrix>,
    relation: Relation,
    name: &str,
) -> Result<(), SymplexError> {
    let n = p.num_vars();
    match (a, b) {
        (None, None) => Ok(()),
        (Some(a), Some(b)) => {
            if a.ncols() != n {
                return Err(invalid(
                    "linprog_matrix",
                    format!("A_{name} has {} columns but c has {n} entries", a.ncols()),
                ));
            }
            let bv = matrix_to_vec(b, &format!("b_{name}"))?;
            if bv.len() != a.nrows() {
                return Err(invalid(
                    "linprog_matrix",
                    format!(
                        "A_{name} has {} rows but b_{name} has {} entries",
                        a.nrows(),
                        bv.len()
                    ),
                ));
            }
            let rows = matrix_to_q(a, &format!("A_{name}"))?;
            for (row, rhs) in rows.into_iter().zip(bv) {
                p.add(row, rhs, relation);
            }
            Ok(())
        }
        _ => Err(invalid(
            "linprog_matrix",
            format!("A_{name} and b_{name} must be given together"),
        )),
    }
}

/// Solve an LP given as [`Matrix`] data with numeric-literal entries.
///
/// * `c`: the objective as an `n×1` column (or `1×n` row) vector;
/// * `a_ub`, `b_ub`: `A_ub x ≤ b_ub` (`m₁×n` and `m₁×1`), both or neither;
/// * `a_eq`, `b_eq`: `A_eq x = b_eq`, both or neither;
/// * bounds are the default `x ≥ 0` (use [`LpProblem`] for anything else).
///
/// Entries are constant-folded with `eval()` first, so `1/2 + 1/3` is
/// accepted; a symbol or transcendental constant is rejected.  Constraints
/// are numbered `≤` rows first, then `=` rows.
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if an entry is not a numeric literal,
/// a matrix is missing its partner, or the shapes are inconsistent.
///
/// # Examples
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::{linprog_matrix, Objective, LpStatus, qi};
///
/// let ctx = Context::new();
/// let c = matrix![ctx, [3], [2]];
/// let a = matrix![ctx, [1, 1], [1, 3]];
/// let b = matrix![ctx, [4], [6]];
/// let sol = linprog_matrix(Objective::Maximize, &c, Some(&a), Some(&b), None, None).unwrap();
/// assert_eq!(sol.status, LpStatus::Optimal);
/// assert_eq!(sol.x_ex(&ctx), vec![ctx.int(4), ctx.int(0)]);
/// assert_eq!(sol.objective, Some(qi(12)));
///
/// let x = ctx.symbol("x");
/// let bad = Matrix::new(vec![vec![x, ctx.int(1)]]).unwrap();
/// assert!(linprog_matrix(Objective::Minimize, &c, Some(&bad), Some(&matrix![ctx, [1]]), None, None).is_err());
/// ```
pub fn linprog_matrix(
    objective: Objective,
    c: &Matrix,
    a_ub: Option<&Matrix>,
    b_ub: Option<&Matrix>,
    a_eq: Option<&Matrix>,
    b_eq: Option<&Matrix>,
) -> Result<LpSolution, SymplexError> {
    let cv = matrix_to_vec(c, "c")?;
    let mut p = match objective {
        Objective::Minimize => LpProblem::minimize(cv),
        Objective::Maximize => LpProblem::maximize(cv),
    };
    add_matrix_block(&mut p, a_ub, b_ub, Relation::Le, "ub")?;
    add_matrix_block(&mut p, a_eq, b_eq, Relation::Eq, "eq")?;
    p.solve()
}

// ═══════════════════════════════════════════════════════════════════════════
// Unit tests
// ═══════════════════════════════════════════════════════════════════════════

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn textbook_max() {
        let sol = LpProblem::maximize(vec![qi(3), qi(2)])
            .le(vec![qi(1), qi(1)], qi(4))
            .le(vec![qi(1), qi(3)], qi(6))
            .solve()
            .unwrap();
        assert_eq!(sol.status, LpStatus::Optimal);
        assert_eq!(sol.x, vec![qi(4), qi(0)]);
        assert_eq!(sol.objective, Some(qi(12)));
    }

    #[test]
    fn infeasible_has_certificate() {
        // x + y ≤ 1 and x + y ≥ 2 with x, y ≥ 0.
        let sol = LpProblem::minimize(vec![qi(1), qi(1)])
            .le(vec![qi(1), qi(1)], qi(1))
            .ge(vec![qi(1), qi(1)], qi(2))
            .solve()
            .unwrap();
        assert_eq!(sol.status, LpStatus::Infeasible);
        let y = sol.farkas.unwrap();
        assert!(!y[0].is_negative());
        assert!(!y[1].is_positive());
        // g = Aᵀy = (y0 + y1, y0 + y1) must be ≥ 0 (lower bounds 0) with
        // inf = 0, and yᵀb = y0 + 2 y1 must be < 0.
        let g = &y[0] + &y[1];
        assert!(!g.is_negative());
        assert!((&y[0] + &(&y[1] * &qi(2))).is_negative());
    }

    #[test]
    fn unbounded_detected() {
        let sol = LpProblem::maximize(vec![qi(1), qi(0)])
            .le(vec![qi(0), qi(1)], qi(1))
            .solve()
            .unwrap();
        assert_eq!(sol.status, LpStatus::Unbounded);
    }

    #[test]
    fn free_variable_negative_optimum() {
        // min x  s.t. x ≥ −3 (as a constraint), x free.
        let sol = LpProblem::minimize(vec![qi(1)])
            .ge(vec![qi(1)], qi(-3))
            .free(0)
            .solve()
            .unwrap();
        assert_eq!(sol.status, LpStatus::Optimal);
        assert_eq!(sol.x, vec![qi(-3)]);
    }

    /// White-box: an artificial that stays basic at level zero after
    /// phase 1 is driven out on a *negative* genuine entry, which flips the
    /// sign of the integer tableau's common denominator `d`.  Phase 2 must
    /// then pivot and read values correctly with `d < 0`.
    #[test]
    fn negative_common_denominator_after_driving_out_artificials() {
        // Row 0: −x − y = 0 (rhs 0, so it is not negated by standardise;
        // both genuine entries are negative, so phase 1 never pivots on
        // them).  Row 1: x + z ≤ 5.  Objective: minimise z, which forces a
        // phase-2 pivot (the slack of row 1 enters).
        let p = LpProblem::minimize(vec![qi(0), qi(0), qi(1)])
            .eq(vec![qi(-1), qi(-1), qi(0)], qi(0))
            .le(vec![qi(1), qi(0), qi(1)], qi(5));
        let Standardized::Ready(sf) = standardize(&p) else {
            panic!("bounds are fine");
        };
        // The same white-box walk on both cell types: the hybrid design's
        // promise is that they take the same path.
        fn walk<T: Cell>(sf: &Standard) -> (Vec<Q>, usize) {
            let mut t = Tableau::<T>::new(sf).expect("fits");
            let (m, n) = (t.m, t.n);
            let phase1 = t.phase1_costs();
            t.set_objective(&phase1).expect("fits");
            assert!(matches!(t.run(n + m).ok().unwrap(), Step::Optimal));
            assert!(t.neg_objective().is_zero(), "feasible");
            assert!(t.basis[0] >= n, "artificial of row 0 still basic");
            assert_eq!(t.d.signum(), std::cmp::Ordering::Greater);
            t.drive_out_artificials().expect("fits");
            assert!(t.basis[0] < n, "artificial driven out");
            assert!(
                t.d.is_negative(),
                "pivot on a negative entry flips the common denominator (d = {:?})",
                t.d
            );
            let pivots_before = t.pivots_done;
            let mut phase2 = vec![Q::zero(); n + m];
            phase2[..n].clone_from_slice(&sf.c);
            t.set_objective(&phase2).expect("fits");
            assert!(matches!(t.run(n).ok().unwrap(), Step::Optimal));
            assert!(t.pivots_done > pivots_before, "phase 2 pivoted with d < 0");
            let z = t.solution();
            assert!(z.iter().all(|v| !v.is_negative()), "z = {z:?}");
            (z, t.pivots_done)
        }
        assert_eq!(walk::<i64>(&sf), walk::<BigInt>(&sf));
        // Cross-check the end-to-end driver on the same problem.
        let sol = p.solve().unwrap();
        assert_eq!(sol.status, LpStatus::Optimal);
        assert_eq!(sol.x, vec![qi(0), qi(0), qi(0)]);
        assert_eq!(sol.objective, Some(qi(0)));
        assert_eq!(sol.duals[1], qi(0));
    }

    /// The 256-bit intermediates of the `i128` cells agree with `BigInt`
    /// on random operands spanning the whole range, including the exact
    /// division's fit/no-fit boundary.
    #[test]
    fn i256_intermediates_match_bigint() {
        let mut state: u128 = 0x9E37_79B9_7F4A_7C15_F39C_C060_5CED_C834;
        let mut next = |bits: u32| -> i128 {
            state ^= state << 13;
            state ^= state >> 7;
            state ^= state << 17;
            let mask = if bits >= 128 {
                u128::MAX
            } else {
                (1u128 << bits) - 1
            };
            let v = (state & mask) as i128;
            if state & (1 << 5) != 0 { -v } else { v }
        };
        for _ in 0..3000 {
            let (a, b, c, d) = (next(120), next(120), next(120), next(120));
            let big = |x: i128| BigInt::from(x);
            // Products and comparisons.
            let ab = I256::mul(a, b);
            let cd = I256::mul(c, d);
            assert_eq!(
                ab.cmp(&cd),
                (big(a) * big(b)).cmp(&(big(c) * big(d))),
                "{a} {b} {c} {d}"
            );
            assert_eq!(
                <i128 as Cell>::cmp_products(&a, &b, &c, &d),
                (big(a) * big(b)).cmp(&(big(c) * big(d)))
            );
            // Exact division: v = d₀·k (both ≤ 60 bits, so v fits), then
            // (v·b)/d₀ = k·b exactly; it fits an i128 for a 64-bit b and
            // usually not for a 120-bit one.
            let d0 = next(60);
            if d0 != 0 {
                let k = next(60);
                let v = d0 * k;
                let b2 = if k & 1 == 0 { next(64) } else { b };
                let exp = big(k) * big(b2);
                assert_eq!(
                    <i128 as Cell>::rescale(&v, &b2, &d0),
                    i128::try_from(&exp).ok(),
                    "{v} {b2} {d0}"
                );
                // The general update (v·p − f·pr)/d with f·pr also a multiple of d₀.
                let f = d0 * next(30);
                let pr = next(64);
                let exp2 = big(k) * big(b2) - big(f) / big(d0) * big(pr);
                assert_eq!(
                    <i128 as Cell>::pivot_update(&v, &b2, &f, &pr, &d0),
                    i128::try_from(&exp2).ok(),
                    "{v} {b2} {f} {pr} {d0}"
                );
            }
            let diff = ab.sub(cd);
            let exp = big(a) * big(b) - big(c) * big(d);
            assert_eq!(diff.to_i128(), i128::try_from(&exp).ok(), "{exp}");
            assert_eq!(diff.is_negative(), Signed::is_negative(&exp));
        }
        // Exact quotients that fill the whole i128 range.
        for &(x, d) in &[
            (i128::MAX, 1i128),
            (i128::MIN, 1),
            (i128::MIN, -1),
            (i128::MAX, i128::MAX),
            (1 << 100, 1 << 40),
        ] {
            let prod = I256::mul(x, d);
            let got = prod.div_exact(d);
            let expected =
                i128::try_from(&(BigInt::from(x) * BigInt::from(d) / BigInt::from(d))).ok();
            assert_eq!(got, expected, "{x} · {d} / {d}");
        }
        // A quotient one past the range is reported as not fitting.
        let over = I256::mul(i128::MAX, 4).sub(I256::from_i128(0));
        assert_eq!(over.div_exact(2), None);
        assert_eq!(I256::mul(i128::MAX, 4).div_exact(4), Some(i128::MAX));
    }

    /// White-box: an LP whose tableau entries outgrow `i64` must fall
    /// back to `BigInt` cells and give the same answer the `BigInt` run
    /// gives on its own; small problems must stay in `i64`.
    #[test]
    fn hybrid_arithmetic_falls_back_to_bigint_on_overflow() {
        // Coefficients around 2^40: a single fraction-free pivot multiplies
        // two of them, so the i64 attempt overflows immediately.
        let big = |k: i64| qi(k) * qi(1 << 40);
        let p = LpProblem::minimize(vec![qi(1), qi(1), qi(1)])
            .ge(vec![big(3), big(1), big(2)], big(7))
            .ge(vec![big(1), big(5), big(1)], big(11))
            .le(vec![big(2), big(1), big(3)], big(40));
        let Standardized::Ready(sf) = standardize(&p) else {
            panic!("bounds are fine");
        };
        assert!(
            matches!(solve_standard::<i64>(&p, &sf), Err(Halt::Overflow)),
            "the i64 attempt must report overflow, not a wrong answer"
        );
        let sol = p.solve().unwrap();
        let via_big = solve_standard::<BigInt>(&p, &sf).ok().unwrap();
        assert_eq!(sol.status, LpStatus::Optimal);
        assert_eq!(sol.x, via_big.x);
        assert_eq!(sol.objective, via_big.objective);
        assert_eq!(sol.duals, via_big.duals);
        // The optimum is exact: 3x + y + 2z ≥ 7, x + 5y + z ≥ 11 (scaled),
        // minimise x + y + z.
        let obj = sol.objective.unwrap();
        assert!(obj.is_positive());
        // A small LP never leaves the i64 path.
        let small = LpProblem::minimize(vec![qi(1), qi(2)])
            .ge(vec![qi(1), qi(1)], qi(1))
            .le(vec![qi(3), qi(1)], qi(6));
        let Standardized::Ready(sf2) = standardize(&small) else {
            panic!("bounds are fine");
        };
        assert!(solve_standard::<i64>(&small, &sf2).is_ok());
    }

    #[test]
    fn malformed_inputs_are_errors() {
        assert!(LpProblem::minimize(vec![]).solve().is_err());
        assert!(
            LpProblem::minimize(vec![qi(1)])
                .le(vec![qi(1), qi(2)], qi(1))
                .solve()
                .is_err()
        );
        assert!(
            LpProblem::minimize(vec![qi(1)])
                .bounds(3, None, None)
                .solve()
                .is_err()
        );
    }
}