Linear Programming

Why Linear Programming Matters

A factory can choose production quantities while respecting machine hours, material limits, and demand caps. If each constraint and the total contribution are linear in the quantities, the model is a linear program, or LP. The same structure appears in transportation, blending, network flows, and many continuous planning models.

Linear programming optimizes continuous variables. It does not make a quantity integral merely because the input coefficients are integers. Its power comes from the geometry of a polyhedron and from dual certificates that can prove an optimum, not just suggest a good point.

This lesson implements primal tableau simplex with deterministic anti-cycling choices and exact rational arithmetic on small integer models. It also explains how revised simplex performs the same basis logic through linear solves, and how interior-point methods approach optimality using primal–dual equations rather than a vertex walk. The code is deliberately a feasible-start Phase II solver: it accepts nonnegative right-hand sides and does not implement Phase I. General feasibility handling is explained separately rather than mislabeled as an infeasibility test.

Core Concepts and Terminology

Use the maximization convention

$$\max_{\mathbf x}\ \mathbf c^{\mathsf T}\mathbf x\qquad\text{subject to }A\mathbf x\leq\mathbf b,\quad\mathbf x\geq\mathbf0.$$

There are $n$ original decision variables and $m$ inequality rows. Our implementation requires integer data, $\mathbf b\geq\mathbf0$, and explicit nonnegativity of every decision variable. The origin is therefore feasible. Negative coefficients in $A$ or $\mathbf c$ are allowed.

Introduce nonnegative slack variables $\mathbf s=\mathbf b-A\mathbf x$, and collect $\mathbf w=(\mathbf x,\mathbf s)$, $M=[A\ I]$, and $\overline{\mathbf c}=(\mathbf c,\mathbf0)$. The standard-form equations are

$$M\mathbf w=\mathbf b,\qquad\mathbf w\geq\mathbf0,\qquad z=\overline{\mathbf c}^{\mathsf T}\mathbf w.$$

A basis is a choice of $m$ linearly independent columns of $M$. Write the corresponding square matrix as $B$. Set nonbasic variables to zero and solve $\mathbf w_B=B^{-1}\mathbf b$. If these basic values are nonnegative, this gives a basic feasible solution. The initial basis consists of the slack columns, so $B=I$ and $\mathbf w_B=\mathbf b$.

A basic variable can equal zero. Such a solution is degenerate: a pivot may change the basis without moving the decision vector or increasing its objective. Bases and geometric vertices are related but are not interchangeable objects.

Tableau convention. Constraint rows store $D=B^{-1}M$ and right-hand side $\mathbf h=B^{-1}\mathbf b$. The objective row represents

$$z+\sum_{j=1}^{n+m}p_jw_j=z_0,\qquad \mathbf p^{\mathsf T}=\mathbf y^{\mathsf T}M-\overline{\mathbf c}^{\mathsf T},\qquad B^{\mathsf T}\mathbf y=\overline{\mathbf c}_B.$$

Thus $p_j$ is minus the reduced cost $r_j=\overline c_j-M_j^{\mathsf T}\mathbf y$. A negative objective-row coefficient means that increasing that nonbasic variable improves a maximization objective. Initially the bottom row is $(-\mathbf c,\mathbf0\mid0)$. Its right-hand side always stores the current objective, with no sign reversal.

The dual of the original maximization model is

$$\min_{\mathbf y}\ \mathbf b^{\mathsf T}\mathbf y\qquad\text{subject to }A^{\mathsf T}\mathbf y\geq\mathbf c,\quad\mathbf y\geq\mathbf0.$$

For feasible $\mathbf x,\mathbf y$, weak duality gives $\mathbf c^{\mathsf T}\mathbf x\leq\mathbf b^{\mathsf T}\mathbf y$. Equality proves both are optimal. In an optimal tableau, the objective-row coefficients in the original slack columns give this $\mathbf y$.

An unbounded maximization LP has feasible points with arbitrarily large objective. A certificate is a feasible point $\mathbf x$ and a direction $\mathbf d\geq\mathbf0$ with $A\mathbf d\leq\mathbf0$ and $\mathbf c^{\mathsf T}\mathbf d\gt0$. Then $\mathbf x+t\mathbf d$ is feasible and improves without limit for $t\geq0$. An infeasible LP has no feasible point; this is a different condition. Every model accepted by this solver is already feasible, so it never returns INFEASIBLE.

For broader geometry, standard-form transformations, and duality, see Goemans's MIT linear-programming notes.

How Linear Programming Works, Step by Step

The implemented algorithm maintains a feasible basis and applies Bland's rule. Variable IDs put original variables first, then slacks; the code uses zero-based IDs, while the sample output uses names such as x1 and s1.

1. Validate the supported model. Check dimensions, integer coefficients, nonnegative right-hand sides, and the pivot cap. Build constraint rows $[A\ I\mid\mathbf b]$ and the objective row $[-\mathbf c\ \mathbf0\mid0]$. Record the initial basis and point.
2. Test optimality. Scan all variable columns in increasing ID order. If no objective coefficient is negative, return OPTIMAL and read the dual vector from the slack columns.
3. Choose the entering variable. Select the smallest variable ID $e$ with $p_e\lt0$. Basic columns have objective coefficient zero in exact arithmetic.
4. Find the first blocking constraint. Let $\mathbf d=D_e$ be the entering column. Increasing the entering variable by $\theta$ changes basic values to $\mathbf h-\theta\mathbf d$. Only rows with $d_i\gt0$ impose upper limits. Minimize $h_i/d_i$ over those rows. If ratios tie, choose the row whose current basic variable has the smallest ID; row position is not the tie-break.
5. Recognize unboundedness if no row blocks. If all $d_i\leq0$, set the entering component of a standard-form direction to $1$, basic components to $-d_i$, and all other nonbasic components to zero. Return its original-variable part as an exact improving ray.
6. Check the pivot cap. If a pivot is required but the permitted count has been reached, return PIVOT_LIMIT. Optimality and unboundedness checks take precedence because neither requires another pivot. A zero cap can therefore still return a certificate.
7. Pivot. Divide the leaving row by its positive pivot entry. Subtract suitable multiples of that normalized row from every other row, including the objective, so the entering column becomes a unit basis column. Replace the leaving variable in the basis list.
8. Record and repeat. Store the entire tableau, original-variable point, current objective, entering and leaving IDs, and ratio. The row numbered $k$ is the state after $k$ completed pivots; its entering, leaving, and ratio fields describe the pivot into that state. Initial fields are NONE and zero.

Bland's entering rule and leaving-variable tie-break prevent cycling in exact arithmetic. They do not require the largest objective gain or the largest reduced cost. The solver records degenerate pivots even when their ratio is zero.

Revised simplex: the same decisions through solves. A large implementation need not form the whole transformed matrix $B^{-1}M$. Maintain the original sparse matrix and a factorization of $B$. Solve $B\mathbf w_B=\mathbf b$ for the basic values and $B^{\mathsf T}\mathbf y=\overline{\mathbf c}_B$ for dual multipliers. Price candidate columns using $r_j=\overline c_j-M_j^{\mathsf T}\mathbf y$. For the chosen entering column, solve $B\mathbf d=M_e$ and perform the same ratio test. Update the basis factorization, with periodic refactorization. Explicitly computing and repeatedly multiplying a stored inverse is generally a poor numerical implementation strategy. Revised simplex is a representation and linear-algebra improvement, not a new optimality criterion.

Interior-point overview. In the standard-form notation above, let the dual slack be $\mathbf q=M^{\mathsf T}\mathbf y-\overline{\mathbf c}$. Primal–dual optimality requires

$$M\mathbf w=\mathbf b,\qquad M^{\mathsf T}\mathbf y-\mathbf q=\overline{\mathbf c},\qquad \mathbf w,\mathbf q\geq\mathbf0,\qquad w_jq_j=0.$$

A central-path approach replaces complementarity by $w_jq_j=\mu$ with $\mu\gt0$, maintaining positive primal variables and dual slacks where the relevant strict-feasibility assumptions hold. Newton steps approximately solve these coupled equations; step control preserves positivity and the target $\mu$ decreases toward zero. These steps generally move through the relative interior rather than visiting adjacent vertices. Practical methods also handle infeasible starting points, scaling, and certificate detection. This is a conceptual overview, not an additional solver hidden inside the code. Boyd and Vandenberghe's interior-point slides develop the barrier and primal–dual viewpoint.

flowchart TD
    A["Validate small integer model and nonnegative RHS"] --> B["Build slack basis and exact tableau; record origin"]
    B --> C{"Any negative objective coefficient?"}
    C -->|No| Z["Return OPTIMAL and exact dual certificate"]
    C -->|Yes| D["Smallest-ID negative coefficient enters"]
    D --> E["Compute ratios only where entering-column entry is positive"]
    E --> F{"Any eligible leaving row?"}
    F -->|No| Y["Return UNBOUNDED with feasible point and improving ray"]
    F -->|Yes| G["Minimum ratio; tie by smallest basic-variable ID"]
    G --> H{"Pivot cap reached?"}
    H -->|Yes| X["Return PIVOT_LIMIT with current feasible point"]
    H -->|No| I["Normalize pivot row; eliminate column in every other row"]
    I --> J["Replace basis variable; record tableau and objective"]
    J --> C

Worked Example with a Full Trace

Consider

$$\max z=3x_1+2x_2,\qquad x_1+x_2\leq4,\quad x_1\leq2,\quad x_2\leq3,\quad x_1,x_2\geq0.$$

Add slacks $s_1,s_2,s_3$ in that row order. The tableau columns below are $(x_1,x_2,s_1,s_2,s_3,\mathrm{rhs})$, with the objective equation last. All entries are exact integers in this example.

Initial state. Basis $(s_1,s_2,s_3)$ has values $(4,2,3)$; original point is $(0,0)$ and $z=0$:

$$T_0=\begin{pmatrix}1&1&1&0&0&4\\1&0&0&1&0&2\\0&1&0&0&1&3\\-3&-2&0&0&0&0\end{pmatrix}.$$

Pivot 1: enter $x_1$, leave $s_2$. The first negative objective coefficient is $-3$, so Bland's rule selects $x_1$. Its constraint column is $(1,1,0)$. Ratios are $4/1=4$ in the $s_1$ row and $2/1=2$ in the $s_2$ row; the zero entry in the $s_3$ row is ineligible. Thus $\theta=2$ and $s_2$ leaves.

The pivot entry is already $1$. Keep the normalized second row; subtract it from the first row; leave the third row unchanged; and add three copies of it to the objective row:

$$T_1=\begin{pmatrix}0&1&1&-1&0&2\\1&0&0&1&0&2\\0&1&0&0&1&3\\0&-2&0&3&0&6\end{pmatrix}.$$

Basis $(s_1,x_1,s_3)$ now has values $(2,2,3)$, so $(x_1,x_2)=(2,0)$. Objective gain is $r_{x_1}\theta=3\cdot2=6$.

Pivot 2: enter $x_2$, leave $s_1$. The only negative objective coefficient is $-2$. The entering column is $(1,0,1)$, and basic values are $(2,2,3)$. Eligible ratios are $2/1=2$ and $3/1=3$; the middle row is skipped. Thus $\theta=2$ and $s_1$ leaves.

The first row is already normalized. Leave the second row unchanged, subtract the first row from the third, and add two copies of the first row to the objective:

$$T_2=\begin{pmatrix}0&1&1&-1&0&2\\1&0&0&1&0&2\\0&0&-1&1&1&1\\0&0&2&1&0&10\end{pmatrix}.$$

Basis $(x_2,x_1,s_3)$ has values $(2,2,1)$. The original point is $(2,2)$, slacks are $(0,0,1)$, and the objective is $6+2\cdot2=10$. No objective coefficient is negative, so the solver returns OPTIMAL.

Independent certificate. The slack-column objective coefficients give $\mathbf y=(2,1,0)$. They are nonnegative and

$$A^{\mathsf T}\mathbf y=(3,2)=\mathbf c,\qquad \mathbf b^{\mathsf T}\mathbf y=4\cdot2+2\cdot1+3\cdot0=10=\mathbf c^{\mathsf T}(2,2).$$

Primal feasibility, dual feasibility, and equality of objectives prove optimality without trusting the pivot history. The certificate is not just a small change in objective or a claim of convergence.

The programs also run three boundary cases. A cap of $1$ on this bounded model returns PIVOT_LIMIT at $(2,0)$ with current value $6$. For $x_1-x_2\leq1$ and objective $x_1+x_2$, simplex reaches $(1,0)$ and returns UNBOUNDED with ray $(1,1)$: $A\mathbf d=0$ and $\mathbf c^{\mathsf T}\mathbf d=2$. For $\max x_1$ subject to $2x_1\leq3$, the exact optimum is $x_1=3/2$ with dual multiplier $1/2$; integer input does not imply an integer LP solution.

Implementation in Python, C++, and Java

Each panel is a complete, independently runnable exact tableau solver for $\max\mathbf c^{\mathsf T}\mathbf x$ with $A\mathbf x\leq\mathbf b$ and $\mathbf x\geq0$. Its signature is simplex(A, b, c, pivot_cap), with a default cap in Python and C++; Java takes the cap explicitly. It accepts $0\leq m\leq3$ rows, $1\leq n\leq6$ variables, integer coefficients in $[-100,100]$, integer right-hand sides in $[0,100]$, and $0\leq\text{pivot cap}\leq2000$. All matrix rows must have length $n$. Python rejects booleans and nonintegers rather than silently converting them. C++ and Java use integer input types.

The nonnegative-RHS requirement gives an immediate feasible slack basis. A negative RHS is rejected as unsupported input, not diagnosed as infeasible. For example, $-x_1\leq-1$, $x_1\geq0$ is feasible, but needs a different start or Phase I and is outside this implementation. Zero rows and even zero inequality constraints are supported; nonnegativity of the decision variables still applies.

Rationals are normalized with positive denominators. Python uses Fraction and Java uses BigInteger. C++ uses a small rational helper backed by the GCC signed 128-bit integer extension, with checked arithmetic that throws on overflow instead of wrapping. The intentionally tiny row and coefficient limits keep simplex intermediates within that representation: tableau entries are ratios of matrix minors of order at most four, and the stated bounds allow a conservative intermediate bound below $10^{28}\lt2^{127}$. Do not raise those limits in C++ without a fresh bound or an arbitrary-precision rational type. Exact arithmetic here removes epsilon tests; it does not turn this dense teaching implementation into a large-scale LP solver.

Run Python 3 with python main.py. Compile the C++ panel using a GCC target with 128-bit integer support: g++ -std=c++17 main.cpp -o lp, then ./lp. Save Java as Main.java and run javac Main.java followed by java Main. No third-party library is needed. Values print as integers or reduced fractions, so there is no decimal formatting or hidden rounding difference.

The result contains a feasible original-variable point and its current objective for every accepted model. OPTIMAL also contains a dual certificate; UNBOUNDED contains an improving ray; PIVOT_LIMIT contains neither. Empty certificate fields print (). The current objective printed for an unbounded run is the value of its finite returned feasible point, not a finite optimum. Each history row owns a complete tableau and basis snapshot. Repeated variable names in different bases refer to the same fixed variable IDs.

The implementations use the same entering order, exact ratio comparisons, leaving-variable tie-break, terminal-check precedence, row-operation order, and pivot cap. They do not implement artificial variables, revised-simplex factorization updates, an interior-point algorithm, integer programming, presolve, or general variable bounds.

from dataclasses import dataclass
from fractions import Fraction as Q


@dataclass
class State:
    iteration: int
    basis: list
    point: list
    value: Q
    entering: int
    leaving: int
    theta: Q
    tableau: list


@dataclass
class Result:
    status: str
    point: list
    value: Q
    dual: list
    ray: list
    pivots: int
    history: list


def simplex(A, b, c, pivot_cap=100):
    m, n = len(b), len(c)
    if not 0 <= m <= 3 or not 1 <= n <= 6 or len(A) != m:
        raise ValueError("supported sizes: 0..3 rows, 1..6 decision variables")
    if any(len(row) != n for row in A):
        raise ValueError("invalid matrix shape")
    if type(pivot_cap) is not int or not 0 <= pivot_cap <= 2000:
        raise ValueError("invalid pivot cap")
    if any(type(v) is not int or not -100 <= v <= 100 for row in A for v in row):
        raise ValueError("A must contain bounded integers")
    if any(type(v) is not int or not 0 <= v <= 100 for v in b):
        raise ValueError("b must contain integers in [0,100]; Phase I is not implemented")
    if any(type(v) is not int or not -100 <= v <= 100 for v in c):
        raise ValueError("c must contain bounded integers")
    N = n + m
    T = [[Q(0) for _ in range(N + 1)] for _ in range(m + 1)]
    for i in range(m):
        for j in range(n):
            T[i][j] = Q(A[i][j])
        T[i][n + i], T[i][N] = Q(1), Q(b[i])
    for j in range(n):
        T[m][j] = Q(-c[j])
    basis = list(range(n, N))
    history = []

    def point():
        x = [Q(0)] * n
        for i, j in enumerate(basis):
            if j < n:
                x[j] = T[i][N]
        return x

    def record(entering, leaving, theta):
        history.append(State(len(history), basis[:], point(), T[m][N],
                             entering, leaving, theta, [row[:] for row in T]))

    def finish(status, dual=None, ray=None):
        return Result(status, point(), T[m][N], dual or [], ray or [],
                      len(history) - 1, history)

    record(-1, -1, Q(0))
    while True:
        entering = next((j for j in range(N) if T[m][j] < 0), -1)
        if entering == -1:
            return finish("OPTIMAL", [T[m][n + i] for i in range(m)])
        row, theta = -1, Q(0)
        for i in range(m):
            if T[i][entering] > 0:
                ratio = T[i][N] / T[i][entering]
                if row == -1 or ratio < theta or (ratio == theta and basis[i] < basis[row]):
                    row, theta = i, ratio
        if row == -1:
            direction = [Q(0)] * N
            direction[entering] = Q(1)
            for i, j in enumerate(basis):
                direction[j] = -T[i][entering]
            return finish("UNBOUNDED", ray=direction[:n])
        if len(history) - 1 == pivot_cap:
            return finish("PIVOT_LIMIT")
        leaving = basis[row]
        pivot = T[row][entering]
        T[row] = [v / pivot for v in T[row]]
        for i in range(m + 1):
            if i != row:
                factor = T[i][entering]
                for j in range(N + 1):
                    T[i][j] = T[i][j] - factor * T[row][j]
        basis[row] = entering
        record(entering, leaving, theta)


def text_q(q):
    return str(q.numerator) if q.denominator == 1 else f"{q.numerator}/{q.denominator}"


def vector(xs):
    return "(" + ",".join(map(str, xs)) + ")"


def show(label, A, b, c, cap):
    r = simplex(A, b, c, cap)
    n = len(c)

    def name(j):
        return "NONE" if j < 0 else f"x{j + 1}" if j < n else f"s{j - n + 1}"

    print(label)
    for s in r.history:
        print(f"k={s.iteration} basis={vector([name(j) for j in s.basis])} "
              f"x={vector([text_q(v) for v in s.point])} z={text_q(s.value)} "
              f"enter={name(s.entering)} leave={name(s.leaving)} theta={text_q(s.theta)}")
    print(f"status={r.status} current_z={text_q(r.value)} pivots={r.pivots} "
          f"dual={vector([text_q(v) for v in r.dual])} ray={vector([text_q(v) for v in r.ray])}")


def main():
    A, b, c = [[1, 1], [1, 0], [0, 1]], [4, 2, 3], [3, 2]
    show("bounded", A, b, c, 100)
    show("pivot-limit", A, b, c, 1)
    show("unbounded", [[1, -1]], [1], [1, 1], 100)
    show("fractional", [[2]], [3], [1], 100)


if __name__ == "__main__":
    main()

#include <algorithm>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;

__extension__ typedef __int128 I;
I add_i(I a, I b) {
    I r; if (__builtin_add_overflow(a, b, &r)) throw overflow_error("rational overflow");
    return r;
}
I sub_i(I a, I b) {
    I r; if (__builtin_sub_overflow(a, b, &r)) throw overflow_error("rational overflow");
    return r;
}
I mul_i(I a, I b) {
    I r; if (__builtin_mul_overflow(a, b, &r)) throw overflow_error("rational overflow");
    return r;
}
I gcd_i(I a, I b) {
    if (a < 0) a = sub_i(0, a);
    if (b < 0) b = sub_i(0, b);
    while (b != 0) { I r = a % b; a = b; b = r; }
    return a;
}
string integer_text(I x) {
    if (x == 0) return "0";
    bool negative = x < 0;
    if (negative) x = sub_i(0, x);
    string s;
    while (x != 0) { s.push_back(char('0' + int(x % 10))); x /= 10; }
    if (negative) s.push_back('-');
    reverse(s.begin(), s.end());
    return s;
}
struct Q {
    I num, den;
    Q(I a = 0, I b = 1) {
        if (b == 0) throw domain_error("zero denominator");
        if (b < 0) { a = sub_i(0, a); b = sub_i(0, b); }
        I g = gcd_i(a, b);
        num = a / g; den = b / g;
    }
    Q operator-() const { return Q(sub_i(0, num), den); }
    Q operator+(const Q& b) const {
        return Q(add_i(mul_i(num, b.den), mul_i(b.num, den)), mul_i(den, b.den));
    }
    Q operator-(const Q& b) const {
        return Q(sub_i(mul_i(num, b.den), mul_i(b.num, den)), mul_i(den, b.den));
    }
    Q operator*(const Q& b) const { return Q(mul_i(num, b.num), mul_i(den, b.den)); }
    Q operator/(const Q& b) const { return Q(mul_i(num, b.den), mul_i(den, b.num)); }
    bool operator<(const Q& b) const { return mul_i(num, b.den) < mul_i(b.num, den); }
    bool operator>(const Q& b) const { return b < *this; }
    bool operator==(const Q& b) const { return num == b.num && den == b.den; }
    string text() const {
        return integer_text(num) + (den == 1 ? "" : "/" + integer_text(den));
    }
};
using Matrix = vector<vector<Q>>;
struct State {
    int iteration;
    vector<int> basis;
    vector<Q> point;
    Q value;
    int entering, leaving;
    Q theta;
    Matrix tableau;
};
struct Result {
    string status;
    vector<Q> point;
    Q value;
    vector<Q> dual, ray;
    int pivots;
    vector<State> history;
};

Result simplex(const vector<vector<int>>& A, const vector<int>& b, const vector<int>& c,
               int pivot_cap = 100) {
    if (b.size() > 3 || c.empty() || c.size() > 6 || A.size() != b.size())
        throw invalid_argument("supported sizes: 0..3 rows, 1..6 decision variables");
    if (pivot_cap < 0 || pivot_cap > 2000) throw invalid_argument("invalid pivot cap");
    for (const auto& row : A) {
        if (row.size() != c.size()) throw invalid_argument("invalid matrix shape");
        for (int v : row) if (v < -100 || v > 100) throw invalid_argument("invalid A");
    }
    for (int v : b) if (v < 0 || v > 100) throw invalid_argument("b needs nonnegative integers");
    for (int v : c) if (v < -100 || v > 100) throw invalid_argument("invalid c");
    int m = int(b.size()), n = int(c.size()), N = n + m;
    Matrix T(m + 1, vector<Q>(N + 1));
    vector<int> basis(m);
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) T[i][j] = Q(A[i][j]);
        T[i][n + i] = Q(1); T[i][N] = Q(b[i]); basis[i] = n + i;
    }
    for (int j = 0; j < n; ++j) T[m][j] = Q(-c[j]);
    vector<State> history;
    auto point = [&]() {
        vector<Q> x(n);
        for (int i = 0; i < m; ++i) if (basis[i] < n) x[basis[i]] = T[i][N];
        return x;
    };
    auto record = [&](int entering, int leaving, Q theta) {
        history.push_back(State{int(history.size()), basis, point(), T[m][N],
                                entering, leaving, theta, T});
    };
    auto finish = [&](const string& status, vector<Q> dual, vector<Q> ray) {
        return Result{status, point(), T[m][N], dual, ray, int(history.size()) - 1, history};
    };
    record(-1, -1, Q(0));
    while (true) {
        int entering = -1;
        for (int j = 0; j < N; ++j) if (T[m][j] < Q(0)) { entering = j; break; }
        if (entering == -1) {
            vector<Q> dual(m);
            for (int i = 0; i < m; ++i) dual[i] = T[m][n + i];
            return finish("OPTIMAL", dual, {});
        }
        int row = -1;
        Q theta;
        for (int i = 0; i < m; ++i) if (T[i][entering] > Q(0)) {
            Q ratio = T[i][N] / T[i][entering];
            if (row == -1 || ratio < theta || (ratio == theta && basis[i] < basis[row])) {
                row = i; theta = ratio;
            }
        }
        if (row == -1) {
            vector<Q> direction(N);
            direction[entering] = Q(1);
            for (int i = 0; i < m; ++i) direction[basis[i]] = -T[i][entering];
            direction.resize(n);
            return finish("UNBOUNDED", {}, direction);
        }
        if (int(history.size()) - 1 == pivot_cap) return finish("PIVOT_LIMIT", {}, {});
        int leaving = basis[row];
        Q pivot = T[row][entering];
        for (int j = 0; j <= N; ++j) T[row][j] = T[row][j] / pivot;
        for (int i = 0; i <= m; ++i) if (i != row) {
            Q factor = T[i][entering];
            for (int j = 0; j <= N; ++j) T[i][j] = T[i][j] - factor * T[row][j];
        }
        basis[row] = entering;
        record(entering, leaving, theta);
    }
}

string vector_text(const vector<Q>& xs) {
    string s = "(";
    for (size_t j = 0; j < xs.size(); ++j) { if (j) s += ','; s += xs[j].text(); }
    return s + ")";
}
void show(const string& label, const vector<vector<int>>& A, const vector<int>& b,
          const vector<int>& c, int cap) {
    Result r = simplex(A, b, c, cap);
    int n = int(c.size());
    auto name = [&](int j) {
        return j < 0 ? string("NONE") : j < n ? "x" + to_string(j + 1) : "s" + to_string(j - n + 1);
    };
    cout << label << '\n';
    for (const auto& s : r.history) {
        string basis = "(";
        for (size_t i = 0; i < s.basis.size(); ++i) { if (i) basis += ','; basis += name(s.basis[i]); }
        cout << "k=" << s.iteration << " basis=" << basis << ") x=" << vector_text(s.point)
             << " z=" << s.value.text() << " enter=" << name(s.entering)
             << " leave=" << name(s.leaving) << " theta=" << s.theta.text() << '\n';
    }
    cout << "status=" << r.status << " current_z=" << r.value.text() << " pivots=" << r.pivots
         << " dual=" << vector_text(r.dual) << " ray=" << vector_text(r.ray) << '\n';
}
int main() {
    vector<vector<int>> A;
    A.push_back({1, 1}); A.push_back({1, 0}); A.push_back({0, 1});
    vector<int> b = {4, 2, 3}, c = {3, 2};
    show("bounded", A, b, c, 100);
    show("pivot-limit", A, b, c, 1);
    vector<vector<int>> ray_A; ray_A.push_back({1, -1});
    show("unbounded", ray_A, {1}, {1, 1}, 100);
    vector<vector<int>> frac_A; frac_A.push_back({2});
    show("fractional", frac_A, {3}, {1}, 100);
    return 0;
}

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class Main {
    static class Q implements Comparable<Q> {
        final BigInteger num, den;
        Q(int x) { this(BigInteger.valueOf(x), BigInteger.ONE); }
        Q(BigInteger n, BigInteger d) {
            if (d.signum() == 0) throw new ArithmeticException("zero denominator");
            if (d.signum() < 0) { n = n.negate(); d = d.negate(); }
            BigInteger g = n.gcd(d);
            num = n.divide(g); den = d.divide(g);
        }
        Q neg() { return new Q(num.negate(), den); }
        Q sub(Q q) { return new Q(num.multiply(q.den).subtract(q.num.multiply(den)), den.multiply(q.den)); }
        Q mul(Q q) { return new Q(num.multiply(q.num), den.multiply(q.den)); }
        Q div(Q q) { return new Q(num.multiply(q.den), den.multiply(q.num)); }
        public int compareTo(Q q) { return num.multiply(q.den).compareTo(q.num.multiply(den)); }
        public String toString() { return num + (den.equals(BigInteger.ONE) ? "" : "/" + den); }
    }
    static final Q ZERO = new Q(0), ONE = new Q(1);
    static Q[][] copy(Q[][] a) {
        Q[][] b = new Q[a.length][];
        for (int i = 0; i < a.length; ++i) b[i] = a[i].clone();
        return b;
    }
    static class State {
        int iteration;
        int[] basis;
        Q[] point;
        Q value;
        int entering, leaving;
        Q theta;
        Q[][] tableau;
        State(int k, int[] b, Q[] x, Q z, int e, int l, Q t, Q[][] a) {
            iteration = k; basis = b.clone(); point = x.clone(); value = z;
            entering = e; leaving = l; theta = t; tableau = copy(a);
        }
    }
    static class Result {
        String status;
        Q[] point;
        Q value;
        Q[] dual, ray;
        int pivots;
        List<State> history;
        Result(String s, Q[] x, Q z, Q[] y, Q[] d, List<State> rows) {
            status = s; point = x.clone(); value = z; dual = y.clone(); ray = d.clone();
            pivots = rows.size() - 1; history = rows;
        }
    }
    static Q[] point(Q[][] T, int[] basis, int n, int N) {
        Q[] x = new Q[n]; Arrays.fill(x, ZERO);
        for (int i = 0; i < basis.length; ++i) if (basis[i] < n) x[basis[i]] = T[i][N];
        return x;
    }
    static Result simplex(int[][] A, int[] b, int[] c, int pivotCap) {
        int m = b.length, n = c.length;
        if (m > 3 || n < 1 || n > 6 || A.length != m)
            throw new IllegalArgumentException("supported sizes: 0..3 rows, 1..6 variables");
        if (pivotCap < 0 || pivotCap > 2000) throw new IllegalArgumentException("invalid pivot cap");
        for (int[] row : A) {
            if (row.length != n) throw new IllegalArgumentException("invalid matrix shape");
            for (int v : row) if (v < -100 || v > 100) throw new IllegalArgumentException("invalid A");
        }
        for (int v : b) if (v < 0 || v > 100) throw new IllegalArgumentException("b needs nonnegative integers");
        for (int v : c) if (v < -100 || v > 100) throw new IllegalArgumentException("invalid c");
        int N = n + m;
        Q[][] T = new Q[m + 1][N + 1];
        for (Q[] row : T) Arrays.fill(row, ZERO);
        int[] basis = new int[m];
        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < n; ++j) T[i][j] = new Q(A[i][j]);
            T[i][n + i] = ONE; T[i][N] = new Q(b[i]); basis[i] = n + i;
        }
        for (int j = 0; j < n; ++j) T[m][j] = new Q(-c[j]);
        List<State> history = new ArrayList<>();
        history.add(new State(0, basis, point(T, basis, n, N), T[m][N], -1, -1, ZERO, T));
        while (true) {
            int entering = -1;
            for (int j = 0; j < N; ++j) if (T[m][j].compareTo(ZERO) < 0) { entering = j; break; }
            if (entering == -1) {
                Q[] dual = new Q[m];
                for (int i = 0; i < m; ++i) dual[i] = T[m][n + i];
                return new Result("OPTIMAL", point(T, basis, n, N), T[m][N], dual, new Q[0], history);
            }
            int row = -1; Q theta = ZERO;
            for (int i = 0; i < m; ++i) if (T[i][entering].compareTo(ZERO) > 0) {
                Q ratio = T[i][N].div(T[i][entering]);
                if (row == -1 || ratio.compareTo(theta) < 0 ||
                    (ratio.compareTo(theta) == 0 && basis[i] < basis[row])) {
                    row = i; theta = ratio;
                }
            }
            if (row == -1) {
                Q[] direction = new Q[N]; Arrays.fill(direction, ZERO); direction[entering] = ONE;
                for (int i = 0; i < m; ++i) direction[basis[i]] = T[i][entering].neg();
                return new Result("UNBOUNDED", point(T, basis, n, N), T[m][N], new Q[0],
                    Arrays.copyOf(direction, n), history);
            }
            if (history.size() - 1 == pivotCap)
                return new Result("PIVOT_LIMIT", point(T, basis, n, N), T[m][N],
                    new Q[0], new Q[0], history);
            int leaving = basis[row]; Q pivot = T[row][entering];
            for (int j = 0; j <= N; ++j) T[row][j] = T[row][j].div(pivot);
            for (int i = 0; i <= m; ++i) if (i != row) {
                Q factor = T[i][entering];
                for (int j = 0; j <= N; ++j) T[i][j] = T[i][j].sub(factor.mul(T[row][j]));
            }
            basis[row] = entering;
            history.add(new State(history.size(), basis, point(T, basis, n, N),
                T[m][N], entering, leaving, theta, T));
        }
    }
    static String vector(Q[] xs) {
        StringBuilder s = new StringBuilder("(");
        for (int i = 0; i < xs.length; ++i) { if (i > 0) s.append(','); s.append(xs[i]); }
        return s.append(')').toString();
    }
    static String name(int j, int n) {
        return j < 0 ? "NONE" : j < n ? "x" + (j + 1) : "s" + (j - n + 1);
    }
    static void show(String label, int[][] A, int[] b, int[] c, int cap) {
        Result r = simplex(A, b, c, cap); int n = c.length;
        System.out.println(label);
        for (State s : r.history) {
            StringBuilder basis = new StringBuilder("(");
            for (int i = 0; i < s.basis.length; ++i) {
                if (i > 0) basis.append(',');
                basis.append(name(s.basis[i], n));
            }
            System.out.println("k=" + s.iteration + " basis=" + basis + ") x=" + vector(s.point) +
                " z=" + s.value + " enter=" + name(s.entering, n) + " leave=" + name(s.leaving, n) +
                " theta=" + s.theta);
        }
        System.out.println("status=" + r.status + " current_z=" + r.value + " pivots=" + r.pivots +
            " dual=" + vector(r.dual) + " ray=" + vector(r.ray));
    }
    public static void main(String[] args) {
        int[][] A = new int[3][];
        A[0] = new int[]{1, 1}; A[1] = new int[]{1, 0}; A[2] = new int[]{0, 1};
        int[] b = {4, 2, 3}, c = {3, 2};
        show("bounded", A, b, c, 100);
        show("pivot-limit", A, b, c, 1);
        int[][] rayA = new int[1][]; rayA[0] = new int[]{1, -1};
        show("unbounded", rayA, new int[]{1}, new int[]{1, 1}, 100);
        int[][] fracA = new int[1][]; fracA[0] = new int[]{2};
        show("fractional", fracA, new int[]{3}, new int[]{1}, 100);
    }
}

Expected standard output, identical in all three languages:

• bounded
• k=0 basis=(s1,s2,s3) x=(0,0) z=0 enter=NONE leave=NONE theta=0
• k=1 basis=(s1,x1,s3) x=(2,0) z=6 enter=x1 leave=s2 theta=2
• k=2 basis=(x2,x1,s3) x=(2,2) z=10 enter=x2 leave=s1 theta=2
• status=OPTIMAL current_z=10 pivots=2 dual=(2,1,0) ray=()
• pivot-limit
• k=0 basis=(s1,s2,s3) x=(0,0) z=0 enter=NONE leave=NONE theta=0
• k=1 basis=(s1,x1,s3) x=(2,0) z=6 enter=x1 leave=s2 theta=2
• status=PIVOT_LIMIT current_z=6 pivots=1 dual=() ray=()
• unbounded
• k=0 basis=(s1) x=(0,0) z=0 enter=NONE leave=NONE theta=0
• k=1 basis=(x1) x=(1,0) z=1 enter=x1 leave=s1 theta=1
• status=UNBOUNDED current_z=1 pivots=1 dual=() ray=(1,1)
• fractional
• k=0 basis=(s1) x=(0) z=0 enter=NONE leave=NONE theta=0
• k=1 basis=(x1) x=(3/2) z=3/2 enter=x1 leave=s1 theta=3/2
• status=OPTIMAL current_z=3/2 pivots=1 dual=(1/2) ray=()

Complexity, Correctness, and Invariants

Write $N=n+m$ for the total number of original and slack variables, and $P$ for the number of completed pivots. A dense tableau contains $(m+1)(N+1)$ rational entries. Construction, snapshots, and a full pivot each take $O((m+1)(N+1))$ rational work; entering scans take $O(N)$, and ratio scans use $O(m)$ exact comparisons and divisions.

The reference program therefore uses $O((P+1)(m+1)(N+1))$ rational operations and the same order of storage when it retains every tableau. Without history, working storage is $O((m+1)(N+1))$. The best case returns immediately after setup and an optimality or unboundedness scan, at $\Theta((m+1)(N+1))$ setup cost. There is no distribution-free “average pivot count.” For unrestricted dimensions, simplex has exponential worst-case pivot behavior; an anti-cycling rule gives finite termination, not a polynomial-time bound. A cap $K$ gives an upper bound $O((K+1)(m+1)(N+1))$ rational operations, while permitting an unfinished result.

For the deliberately tiny accepted dimensions, there are at most $\binom93=84$ candidate basis sets. Exact Bland pivots do not revisit a basis, so no more than $83$ pivots are needed before classification; the default cap of $100$ exceeds that bound. This small-model fact is not a general LP complexity result.

Rational field operations are not unit-cost operations on arbitrary-size input. Numerator and denominator growth, integer multiplication, division, and greatest-common-divisor reduction matter. Python and Java use arbitrary-precision arithmetic; the C++ helper uses the bounded, checked representation described above. Quoting only a count of tableau cells omits bit complexity.

Feasibility and equation invariant. Initially the tableau is the original slack system and the basic values are nonnegative. Suppose the current tableau represents the same equations with an identity basis, nonnegative basic RHS, and objective $z+\mathbf p^{\mathsf T}\mathbf w=z_0$. Increasing an entering nonbasic variable to $\theta$ makes basic values $h_i-\theta d_i$. The minimum eligible ratio is precisely the largest allowed step before some positive-direction row hits zero. All other basic values remain nonnegative.

Normalizing and eliminating are invertible row operations, so they preserve the feasible equality system and objective equation. The entering column becomes a unit column, the updated basis has the new feasible values, and the original-variable point extracted from it is feasible. Objective change is $(-p_e)\theta\geq0$. A zero ratio is permitted: feasibility is preserved even though the objective does not increase. Bland's rule supplies the separate anti-cycling argument.

Optimality certificate. When every $p_j\geq0$, the current objective equation bounds $z\leq z_0$ for every nonnegative feasible $\mathbf w$. Equivalently, $\mathbf p=M^{\mathsf T}\mathbf y-\overline{\mathbf c}\geq0$. Slack columns imply $\mathbf y\geq0$, and original columns imply $A^{\mathsf T}\mathbf y\geq\mathbf c$. The current basic point has $\mathbf c^{\mathsf T}\mathbf x=z_0=\mathbf b^{\mathsf T}\mathbf y$. Weak duality proves it is globally optimal.

Unboundedness certificate. With $p_e\lt0$ and no positive direction entry, the constructed standard-form direction has nonnegative components and satisfies $M\mathbf d_{\rm full}=0$. Its original-variable part obeys $\mathbf d\geq0$ and $A\mathbf d\leq0$, because the slack direction is also nonnegative. Moreover $\mathbf c^{\mathsf T}\mathbf d=-p_e\gt0$. The returned feasible point plus any nonnegative multiple of that direction is feasible and improves without limit.

These are exact certificates, not convergence heuristics. PIVOT_LIMIT establishes neither optimality nor unboundedness. The exact program also makes a meaningful distinction between a repeated objective value and a repeated basis.

Revised simplex cost. Dense basis solves take $O(m^2)$ after a factorization; a fresh dense factorization costs $O(m^3)$. Pricing all columns typically adds $O(mN)$ arithmetic, with sparse structure and pricing policies changing the practical cost. Sparse revised simplex retains the original matrix and basis factors rather than a full dense transformed tableau. Factor fill, stability, update strategy, and periodic refactorization are central engineering issues.

Interior-point cost and termination. Suitable path-following LP algorithms have polynomial complexity guarantees under their stated input and initialization assumptions. A practical primal–dual iteration is often dominated by a sparse linear-system factorization, not by scanning every neighboring vertex. In the standard-form central-path equations, the duality gap is

$$\mathbf b^{\mathsf T}\mathbf y-\overline{\mathbf c}^{\mathsf T}\mathbf w=\mathbf w^{\mathsf T}\mathbf q=N\mu$$

when primal and dual equations are satisfied and $w_jq_j=\mu$ for every component. Thus decreasing complementarity targets has a direct optimality meaning. Approximate implementations check primal residuals, dual residuals, and gap together; a small gap alone with infeasible vectors is insufficient. General infeasibility and unboundedness handling requires additional machinery, such as appropriate homogeneous formulations. An interior-point iterate need not be a vertex; crossover can construct a basis when that is desired.

Variants, Pitfalls, and Optimizations

Feasible start versus general feasibility. A negative RHS does not mean the LP is infeasible. A general two-phase method first constructs an auxiliary problem with artificial variables, minimizing their sum. A positive auxiliary optimum proves the original equalities cannot be satisfied with the required nonnegativity. A zero optimum permits Phase II after artificial variables are removed, including pivoting out zero-valued artificial basics or deleting redundant rows when appropriate. Restore and canonicalize the original objective before continuing. Multiplying a constraint by $-1$ reverses its inequality direction; adding a slack with the wrong sign changes the problem. Cornell's LP course, Lectures 6–10 covers this transition and general forms.

Respect the tableau sign convention. This maximization solver enters a column with a negative bottom-row coefficient because the row stores minus reduced costs. A different book can use the opposite sign convention. Derive the objective equation before copying a pivot rule. The RHS here is the actual current objective, not its negative.

Use the correct ratio test and tie-break. Ratios are considered only for strictly positive entering-column entries, and the minimum preserves feasibility. Zero ratios are valid degenerate pivots, not a reason to divide by zero or discard the row. Ties are broken by the basic variable's fixed ID, which can differ from the row number after pivots. Choosing the largest ratio can move to a point with a larger objective that violates the constraints.

Degeneracy, cycling, and multiple optima are different. A zero basic variable signals degeneracy. Cycling concerns repeated bases under a pivot rule. Multiple optima concern distinct feasible points with the same optimal value. At an optimal tableau, a zero reduced-cost nonbasic column can indicate an alternative optimal movement, but degeneracy may make a particular zero-cost pivot fail to produce a different point. Do not declare termination merely because one pivot leaves the objective unchanged. Exact Bland pivots prevent cycling here; numerical solvers also need tolerance and stability policies.

LP relaxation is not integer programming. The fractional sample has integer inputs and an optimum $3/2$. If a production variable must be an integer, the LP is only a relaxation unless special structure supplies integrality. Rounding a fractional solution can violate constraints or lose optimality. Network incidence structures can have integrality properties under suitable hypotheses; arbitrary integer matrices do not.

Numerical and modeling practice. Large solvers ordinarily use floating-point sparse factorizations, scaling, presolve, and refactorization. Near-zero reduced costs or pivot entries then need numerical policies, and reported primal/dual residuals must be checked. Exact rational arithmetic is valuable for small demonstrations and independent verification, while bit growth can make it expensive at scale. Polynomial bit complexity and strongly polynomial arithmetic-operation bounds are different notions; coefficient bit length matters in the former.

Variable lower bounds can be shifted, free variables split into a difference of nonnegative variables, and equalities handled in standard form. These transformations change dimensions and may change the feasible starting basis. The reference solver does not silently perform them. Revised simplex warm starts can reuse a basis after suitable model changes, while interior-point solutions may require crossover to obtain a useful basic representation.

Practice Problems and Real-World Use

Start by reproducing the bounded primal and its dual certificate, then verify the unbounded ray without using the tableau. These checks separate optimization claims from implementation details. For larger models, formulate units and bounds carefully and inspect solver status, primal/dual residuals, and certificate information before using the result.

The exercises below connect vertex geometry, exact linear algebra, and the practical representations introduced in the lesson. Keep continuous and integer versions of a model distinct.