diff --git a/python/cuopt/cuopt/linear_programming/problem.py b/python/cuopt/cuopt/linear_programming/problem.py
index 49885acf52..cae15282c4 100644
--- a/python/cuopt/cuopt/linear_programming/problem.py
+++ b/python/cuopt/cuopt/linear_programming/problem.py
@@ -48,6 +48,11 @@ class CType(str, Enum):
     GE = "G"
     EQ = "E"
 
+    @property
+    def symbol(self):
+        """Algebraic symbol used when printing constraints."""
+        return {CType.LE: "<=", CType.GE: ">=", CType.EQ: "=="}[self]
+
 
 LE = CType.LE
 GE = CType.GE
@@ -335,6 +340,131 @@ def __eq__(self, other):
             case _:
                 raise ValueError("Unsupported operation")
 
+    def _display_name(self):
+        if self.VariableName:
+            return self.VariableName
+        if self.index >= 0:
+            # Same name an unnamed variable gets on export.
+            return f"C{self.index}"
+        return ""
+
+    def __str__(self):
+        return self._display_name() or repr(self)
+
+    def __repr__(self):
+        vtype = self.VariableType
+        if isinstance(vtype, (bytes, bytearray)):
+            # VariableType is not normalized, see #1736.
+            vtype = vtype.decode()
+        return (
+            f"<cuopt.Variable {self._display_name()!r} (index={self.index}), "
+            f"type={VType(vtype).name}, bounds=[{self.LB}, {self.UB}], "
+            f"value={self.Value}>"
+        )
+
+
+# Maximum number of terms rendered when stringifying a linear or quadratic
+# expression. Beyond this, the head is shown followed by a ``... (N more
+# terms)`` marker so that printing a model with thousands of terms stays
+# readable in a REPL or notebook instead of flooding the output. Set to
+# ``None`` to disable truncation entirely.
+_MAX_DISPLAY_TERMS = 10
+
+
+def _same_variable(var1, var2):
+    """Identity/index comparison; Variable.__eq__ builds a Constraint."""
+    return var1 is var2 or (var1.index >= 0 and var1.index == var2.index)
+
+
+class _ExprBuilder:
+    """Build an algebraic string from a sequence of terms.
+
+    The first term is emitted without a sign; subsequent terms are joined
+    with ' + ' or ' - ' separators. A coefficient of 1.0 or -1.0 is
+    elided, so '1.0 * x' becomes 'x' and '-1.0 * x' becomes '-x'.
+
+    When ``max_terms`` is set, only the first ``max_terms`` non-zero terms
+    are rendered; any remaining terms are counted and summarized as a
+    trailing ``... (N more terms)`` marker. This keeps the output bounded
+    for expressions with very many terms. ``max_terms=None`` (the default)
+    renders every term.
+    """
+
+    def __init__(self, max_terms=None):
+        self.parts = []
+        self.max_terms = max_terms
+        # Non-zero terms seen so far (rendered + hidden).
+        self.n_terms = 0
+        # Non-zero terms omitted because the cap was reached.
+        self.n_hidden = 0
+
+    def add_linear(self, coef, var):
+        """Add a linear term ``coef * var``."""
+        if coef == 0.0:
+            return
+        var_str = str(var)
+        if coef == 1.0:
+            self._append(var_str, negative=False)
+        elif coef == -1.0:
+            self._append(var_str, negative=True)
+        else:
+            self._append(f"{abs(coef)} * {var_str}", negative=coef < 0)
+
+    def add_quadratic(self, coef, var1, var2):
+        """Add a quadratic term ``coef * var1 * var2``."""
+        if coef == 0.0:
+            return
+        v1_str = str(var1)
+        v2_str = str(var2)
+        if _same_variable(var1, var2):
+            term_str = f"{v1_str}^2"
+        elif v1_str <= v2_str:
+            term_str = f"{v1_str} * {v2_str}"
+        else:
+            term_str = f"{v2_str} * {v1_str}"
+        if coef == 1.0:
+            self._append(term_str, negative=False)
+        elif coef == -1.0:
+            self._append(term_str, negative=True)
+        else:
+            self._append(f"{abs(coef)} * {term_str}", negative=coef < 0)
+
+    def add_constant(self, value):
+        """Add a constant term."""
+        if value == 0.0:
+            return
+        self._append(f"{abs(value)}", negative=value < 0)
+
+    def _append(self, term, negative):
+        self.n_terms += 1
+        if self.max_terms is not None and self.n_terms > self.max_terms:
+            # Past the cap: count the term but don't render it.
+            self.n_hidden += 1
+            return
+        if not self.parts:
+            self.parts.append(f"-{term}" if negative else term)
+        else:
+            self.parts.append(f" - {term}" if negative else f" + {term}")
+
+    def build(self):
+        if not self.parts and not self.n_hidden:
+            return "0.0"
+        result = "".join(self.parts)
+        if self.n_hidden:
+            plural = "term" if self.n_hidden == 1 else "terms"
+            marker = f"... ({self.n_hidden} more {plural})"
+            result = f"{result} + {marker}" if result else marker
+        return result
+
+
+def _format_linear(vars, coeffs, constant, max_terms=None):
+    """Format a linear expression as an algebraic string."""
+    builder = _ExprBuilder(max_terms=max_terms)
+    for var, coef in zip(vars, coeffs):
+        builder.add_linear(coef, var)
+    builder.add_constant(constant)
+    return builder.build()
+
 
 class QuadraticExpression:
     """
@@ -889,6 +1019,25 @@ def __ge__(self, other):
     def __eq__(self, other):
         raise ValueError("Equality constraints are not supported.")
 
+    def __str__(self):
+        builder = _ExprBuilder(max_terms=_MAX_DISPLAY_TERMS)
+        if self.qmatrix is not None:
+            for row, col, val in zip(
+                self.qmatrix.row, self.qmatrix.col, self.qmatrix.data
+            ):
+                if val == 0.0:
+                    continue
+                builder.add_quadratic(val, self.qvars[row], self.qvars[col])
+        for v1, v2, coef in zip(self.qvars1, self.qvars2, self.qcoefficients):
+            builder.add_quadratic(coef, v1, v2)
+        for var, coef in zip(self.vars, self.coefficients):
+            builder.add_linear(coef, var)
+        builder.add_constant(self.constant)
+        return builder.build()
+
+    def __repr__(self):
+        return f"<cuopt.QuadraticExpression: {self}>"
+
 
 def _quadratic_expression_to_qcmatrix(expr, rhs):
     """Build QCMATRIX COO data for a quadratic row ``expr`` sense ``rhs``.
@@ -1280,6 +1429,17 @@ def __eq__(self, other):
                 expr = self - other
                 return Constraint(expr, EQ, 0.0)
 
+    def __str__(self):
+        return _format_linear(
+            self.vars,
+            self.coefficients,
+            self.constant,
+            max_terms=_MAX_DISPLAY_TERMS,
+        )
+
+    def __repr__(self):
+        return f"<cuopt.LinearExpression: {self}>"
+
 
 class Constraint:
     """
@@ -1341,7 +1501,12 @@ def __init__(self, expr, sense, rhs, name=""):
             self.rhs_value = rhs_value
             self.RHS = rhs_value
             self.vindex_coeff_dict = {}
-            self.vars = expr.vars
+            # expr.vars holds only the linear terms; id() because Variable
+            # overrides __eq__ and is unhashable.
+            seen = {}
+            for var in (*expr.vars, *expr.qvars1, *expr.qvars2, *expr.qvars):
+                seen.setdefault(id(var), var)
+            self.vars = list(seen.values())
             return
 
         self.is_quadratic = False
@@ -1397,6 +1562,28 @@ def compute_slack(self):
 
         return self.RHS - lhs
 
+    def __str__(self):
+        # Rendered from the data the constraint stores for the solver, so
+        # the output is normalized: duplicate terms are merged and any
+        # expression constant is folded into the right-hand side.
+        builder = _ExprBuilder(max_terms=_MAX_DISPLAY_TERMS)
+        index_to_var = {v.index: v for v in self.vars}
+        if self.is_quadratic:
+            for row, col, val in zip(self.rows, self.cols, self.vals):
+                builder.add_quadratic(
+                    val, index_to_var[row], index_to_var[col]
+                )
+            for idx, val in zip(self.linear_indices, self.linear_values):
+                builder.add_linear(val, index_to_var[idx])
+        else:
+            for idx, coeff in self.vindex_coeff_dict.items():
+                builder.add_linear(coeff, index_to_var[idx])
+        return f"{builder.build()} {CType(self.Sense).symbol} {self.RHS}"
+
+    def __repr__(self):
+        name = self.ConstraintName if self.ConstraintName else "<unnamed>"
+        return f"<cuopt.Constraint '{name}': {self}>"
+
 
 class Problem:
     """
@@ -1760,9 +1947,15 @@ def updateConstraint(self, constr, coeffs=None, rhs=None):
                 )
             if isinstance(coeffs, dict):
                 coeffs = coeffs.items()
+            new_vars = []
             for var, coeff in coeffs:
                 idx = var.index
+                if idx not in constr.vindex_coeff_dict:
+                    new_vars.append(var)
                 constr.vindex_coeff_dict[idx] = coeff
+            if new_vars:
+                # constr.vars aliases the expression's list; rebind it.
+                constr.vars = constr.vars + new_vars
             if rhs is not None:
                 constr.RHS = rhs
         else:
@@ -2223,3 +2416,54 @@ def solve(self, settings=solver_settings.SolverSettings()):
         # Post Solve
         self.populate_solution(solution)
         return solution
+
+    def __repr__(self):
+        name = self.Name if self.Name else "<unnamed>"
+        return (
+            f"<cuopt.Problem '{name}' "
+            f"({len(self.vars)} vars, {len(self.constrs)} constrs, "
+            f"IsMIP={self.IsMIP})>"
+        )
+
+    def __str__(self):
+        lines = []
+        name = self.Name if self.Name else "<unnamed>"
+        lines.append(f"Problem: {name}")
+        sense_str = "MINIMIZE" if self.ObjSense == MINIMIZE else "MAXIMIZE"
+        lines.append(f"  Objective: {sense_str}")
+
+        n_cont = 0
+        n_int = 0
+        n_semi = 0
+        for v in self.vars:
+            t = v.VariableType
+            if isinstance(t, (bytes, bytearray)):
+                t = t.decode()
+            if t in ("I", VType.INTEGER):
+                n_int += 1
+            elif t in ("S", VType.SEMI_CONTINUOUS):
+                n_semi += 1
+            else:
+                n_cont += 1
+        lines.append(
+            f"  Variables: {len(self.vars)} "
+            f"(continuous={n_cont}, integer={n_int}, "
+            f"semi-continuous={n_semi})"
+        )
+
+        n_linear = sum(1 for c in self.constrs if not c.is_quadratic)
+        n_quad = sum(1 for c in self.constrs if c.is_quadratic)
+        lines.append(
+            f"  Constraints: {len(self.constrs)} "
+            f"(linear={n_linear}, quadratic={n_quad})"
+        )
+        lines.append(f"  Non-zeros: {self.NumNZs}")
+
+        if self.solved:
+            status = self.Status
+            if hasattr(status, "name"):
+                status = status.name
+            lines.append(f"  Status: {status}")
+            lines.append(f"  Objective value: {self.ObjValue}")
+
+        return "\n".join(lines)
diff --git a/python/cuopt/cuopt/tests/linear_programming/test_python_API.py b/python/cuopt/cuopt/tests/linear_programming/test_python_API.py
index 964809553f..f701d094dc 100644
--- a/python/cuopt/cuopt/tests/linear_programming/test_python_API.py
+++ b/python/cuopt/cuopt/tests/linear_programming/test_python_API.py
@@ -18,7 +18,9 @@
     MINIMIZE,
     SEMI_CONTINUOUS,
     CType,
+    LinearExpression,
     Problem,
+    Variable,
     VType,
     sense,
     QuadraticExpression,
@@ -160,6 +162,43 @@ def test_constraint_duplicate_terms_slack():
     assert c.compute_slack() == pytest.approx(6.0)
 
 
+def test_updateConstraint_tracks_new_variables():
+    """Variables added via updateConstraint end up in Constraint.vars."""
+    prob = Problem()
+    x1 = prob.addVariable(name="x1")
+    x2 = prob.addVariable(name="x2")
+    c = prob.addConstraint(2 * x1 <= 10)
+    assert [v.index for v in c.vars] == [0]
+
+    prob.updateConstraint(c, coeffs=[(x2, 4.0)])
+    assert [v.index for v in c.vars] == [0, 1]
+    x1.Value = 1.0
+    x2.Value = 2.0
+    assert c.compute_slack() == pytest.approx(0.0)
+
+
+def test_updateConstraint_does_not_mutate_expression():
+    """The expression a constraint was built from is left unchanged."""
+    prob = Problem()
+    a = prob.addVariable(name="a")
+    b = prob.addVariable(name="b")
+    expr = 2 * a
+    c = prob.addConstraint(expr <= 10)
+    prob.updateConstraint(c, coeffs=[(b, 1.0)])
+    assert len(expr.vars) == 1
+    assert len(expr.coefficients) == 1
+
+
+def test_constraint_vars_includes_quadratic_only_variables():
+    """Variables that appear only in quadratic terms are in Constraint.vars."""
+    prob = Problem()
+    x = prob.addVariable(name="x")
+    y = prob.addVariable(name="y")
+    c = prob.addConstraint(x * x + 2 * x * y <= 4)
+    assert c.is_quadratic
+    assert [v.index for v in c.vars] == [0, 1]
+
+
 def test_semi_continuous_variable():
     prob = Problem("Semi-continuous")
     x = prob.addVariable(lb=5.0, ub=10.0, vtype=SEMI_CONTINUOUS, name="x")
@@ -851,3 +890,193 @@ def test_quadratic_matrix_2():
     assert x2.getValue() == pytest.approx(0.0000000, abs=1e-3)
     assert x3.getValue() == pytest.approx(0.1092896, abs=1e-3)
     assert problem.ObjValue == pytest.approx(3.715847, abs=1e-3)
+
+
+def test_str_and_repr():
+    """Verify algebraic __str__ and detailed __repr__ for LP API classes."""
+    prob = Problem("str_repr_test")
+
+    # === Variable ===
+    x = prob.addVariable(lb=0.0, ub=10.0, vtype=VType.CONTINUOUS, name="x")
+    y = prob.addVariable(lb=0.0, ub=5.0, vtype=VType.INTEGER, name="y")
+    z = prob.addVariable()
+
+    # __str__: with name returns the name
+    assert str(x) == "x"
+    assert str(y) == "y"
+    # __str__: without name falls back to C{index}
+    assert str(z) == "C2"
+    # __str__: outside a problem there is neither name nor index
+    assert str(Variable()).startswith("<cuopt.Variable '' (index=-1), ")
+
+    # __repr__: detailed summary
+    r = repr(x)
+    assert "cuopt.Variable" in r
+    assert "'x'" in r
+    assert "index=0" in r
+    assert "type=CONTINUOUS" in r
+    assert "bounds=[0.0, 10.0]" in r
+    assert "value=nan" in r
+
+    r = repr(y)
+    assert "type=INTEGER" in r
+    assert "bounds=[0.0, 5.0]" in r
+
+    r = repr(z)
+    assert "'C2'" in r
+    assert "index=2" in r
+
+    # === LinearExpression ===
+    expr1 = 2 * x + 3 * y
+    expr2 = expr1 - 5
+    expr3 = -x + 2.5
+
+    # __str__
+    assert str(expr1) == "2.0 * x + 3.0 * y"
+    assert str(expr2) == "2.0 * x + 3.0 * y - 5.0"
+    assert str(expr3) == "-x + 2.5"
+    # Empty expression collapses to 0.0
+    assert str(LinearExpression([], [], 0.0)) == "0.0"
+    # Constant-only expression
+    assert str(LinearExpression([], [], 3.0)) == "3.0"
+    assert str(LinearExpression([], [], -3.0)) == "-3.0"
+
+    # __repr__
+    assert repr(expr1) == "<cuopt.LinearExpression: 2.0 * x + 3.0 * y>"
+
+    # === QuadraticExpression ===
+    qexpr1 = x * x
+    qexpr2 = qexpr1 + 2 * x * y + 3 * x
+    qexpr3 = -x * x + 0.5 * y * y + x * y
+
+    # __str__
+    assert str(qexpr1) == "x^2"
+    assert str(qexpr2) == "x^2 + 2.0 * x * y + 3.0 * x"
+    assert str(qexpr3) == "-x^2 + 0.5 * y^2 + x * y"
+    # Empty quadratic expression
+    assert str(QuadraticExpression()) == "0.0"
+
+    # __str__: two distinct variables sharing a name are not a square
+    dup_prob = Problem()
+    x_a = dup_prob.addVariable(name="x")
+    x_b = dup_prob.addVariable(name="x")
+    assert str(x_a * x_b) == "x * x"
+    assert str(x_a * x_a) == "x^2"
+
+    # __repr__
+    assert repr(qexpr1) == "<cuopt.QuadraticExpression: x^2>"
+
+    # === Constraint ===
+    c1 = 2 * x + 3 * y <= 10
+    c2 = x - y >= 0
+    c3 = x + 1 == 5
+    prob.addConstraint(c1, name="c1")
+    prob.addConstraint(c2, name="c2")
+    prob.addConstraint(c3, name="c3")
+
+    # __str__: shows the normalized form the solver holds (duplicate terms
+    # merged, expression constants folded into the right-hand side)
+    assert str(c1) == "2.0 * x + 3.0 * y <= 10.0"
+    assert str(c2) == "x - y >= 0.0"
+    assert str(c3) == "x == 4.0"
+
+    # __str__: unnamed constraint
+    c_anon = 2 * x + 3 * y <= 10
+    assert "2.0 * x + 3.0 * y <= 10.0" in str(c_anon)
+
+    # __str__: duplicate terms are merged
+    c_dup = 2 * x + 3 * x <= 5
+    assert str(c_dup) == "5.0 * x <= 5.0"
+
+    # __str__: quadratic constraint rendered from its QCMATRIX data
+    c_quad = x * x + 2 * x * y <= 4
+    assert str(c_quad) == "x^2 + 2.0 * x * y <= 4.0"
+
+    # __str__: reflects updateConstraint (no stale expression data)
+    prob_upd = Problem("upd_test")
+    a = prob_upd.addVariable(name="a")
+    b = prob_upd.addVariable(name="b")
+    c_upd = prob_upd.addConstraint(2 * a + 3 * b <= 10, name="c_upd")
+    prob_upd.updateConstraint(c_upd, coeffs=[(a, 7.0)], rhs=20.0)
+    assert str(c_upd) == "7.0 * a + 3.0 * b <= 20.0"
+    # __str__: variable introduced by updateConstraint
+    c_new = prob_upd.addConstraint(2 * a <= 10, name="c_new")
+    prob_upd.updateConstraint(c_new, coeffs=[(b, 4.0)])
+    assert str(c_new) == "2.0 * a + 4.0 * b <= 10.0"
+
+    # __repr__
+    assert repr(c1) == "<cuopt.Constraint 'c1': 2.0 * x + 3.0 * y <= 10.0>"
+    assert repr(c2) == "<cuopt.Constraint 'c2': x - y >= 0.0>"
+    assert repr(c3) == "<cuopt.Constraint 'c3': x == 4.0>"
+
+    # === Problem ===
+    # __repr__
+    r = repr(prob)
+    assert "cuopt.Problem" in r
+    assert "str_repr_test" in r
+    assert "3 vars" in r
+    assert "3 constrs" in r
+    assert "IsMIP=True" in r  # y is integer
+
+    # __str__: before solve
+    s = str(prob)
+    assert "str_repr_test" in s
+    assert "MINIMIZE" in s
+    assert "Variables: 3" in s
+    assert "continuous=2" in s
+    assert "integer=1" in s
+    assert "semi-continuous=0" in s
+    assert "Constraints: 3" in s
+    assert "linear=3" in s
+    assert "quadratic=0" in s
+    assert "Non-zeros: 5" in s
+    # No status before solve
+    assert "Status:" not in s
+    assert "Objective value:" not in s
+
+    # Unnamed problem
+    empty_prob = Problem()
+    assert "Problem: <unnamed>" in str(empty_prob)
+    assert "<cuopt.Problem '<unnamed>'" in repr(empty_prob)
+
+
+def test_problem_str_after_solve():
+    """Problem.__str__ reports status and objective value once solved."""
+    prob = Problem("solved")
+    x = prob.addVariable(lb=0.0, ub=1.0, name="x")
+    prob.setObjective(x, sense=MINIMIZE)
+    settings = SolverSettings()
+    settings.set_parameter("time_limit", 5)
+    prob.solve(settings)
+    s = str(prob)
+    assert "Status: Optimal" in s
+    assert f"Objective value: {prob.ObjValue}" in s
+
+
+def test_str_truncation_large_expression():
+    """Expressions past the display cap end in a "more terms" marker."""
+    prob = Problem()
+    xs = [prob.addVariable(name=f"x{i}") for i in range(12)]
+
+    expr = xs[0] + xs[1]
+    for x in xs[2:]:
+        expr = expr + x
+    assert str(expr) == (
+        "x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + ... (2 more terms)"
+    )
+    assert repr(expr) == f"<cuopt.LinearExpression: {expr}>"
+
+    at_cap = xs[0] + xs[1]
+    for x in xs[2:10]:
+        at_cap = at_cap + x
+    assert str(at_cap) == "x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9"
+    assert str(at_cap + xs[10]).endswith("+ x9 + ... (1 more term)")
+
+    qexpr = xs[0] * xs[0]
+    for x in xs[1:]:
+        qexpr = qexpr + x * x
+    assert str(qexpr) == (
+        "x0^2 + x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2 + x9^2"
+        " + ... (2 more terms)"
+    )
+    assert repr(qexpr) == f"<cuopt.QuadraticExpression: {qexpr}>"