diff --git a/libsolutil/BooleanLP.cpp b/libsolutil/BooleanLP.cpp
index 3fd5a22ea..a5db6f658 100644
--- a/libsolutil/BooleanLP.cpp
+++ b/libsolutil/BooleanLP.cpp
@@ -136,7 +136,7 @@ void BooleanLPSolver::addAssertion(Expression const& _expr)
 			{
 				LinearExpression data = *left - *right;
 				data[0] *= -1;
-				Constraint c{move(data), _expr.name == "=", {}};
+				Constraint c{move(data), _expr.name == "="};
 				if (!tryAddDirectBounds(c))
 					state().fixedConstraints.emplace_back(move(c));
 				cout << "Added as fixed constraint" << endl;
@@ -186,7 +186,7 @@ void BooleanLPSolver::addAssertion(Expression const& _expr)
 		{
 			LinearExpression data = *left - *right;
 			data[0] *= -1;
-			Constraint c{move(data), _expr.name == "=", {}};
+			Constraint c{move(data), _expr.name == "="};
 			if (!tryAddDirectBounds(c))
 				state().fixedConstraints.emplace_back(move(c));
 		}
@@ -220,10 +220,24 @@ pair<CheckResult, vector<string>> BooleanLPSolver::check(vector<Expression> cons
 
 	std::vector<std::string> booleanVariables;
 	std::vector<Clause> clauses = state().clauses;
-	SolvingState lpState;
+
+	// TODO we start building up a new set of solver
+	// for each query, but we should also keep some
+	// kind of cache across queries.
+	std::vector<std::pair<size_t, LPSolver>> lpSolvers;
+	lpSolvers.emplace_back(0, LPSolver{});
+	LPSolver& lpSolver = lpSolvers.back().second;
+
 	for (auto&& [index, bound]: state().bounds)
-		resizeAndSet(lpState.bounds, index, bound);
-	lpState.constraints = state().fixedConstraints;
+	{
+		if (bound.lower)
+			lpSolver.addLowerBound(index, *bound.lower);
+		if (bound.upper)
+			lpSolver.addUpperBound(index, *bound.upper);
+	}
+	for (Constraint const& c: state().fixedConstraints)
+		lpSolver.addConstraint(c);
+
 	// TODO this way, it will result in a lot of gaps in both sets of variables.
 	// should we compress them and store a mapping?
 	// Is it even a problem if the indices overlap?
@@ -231,23 +245,9 @@ pair<CheckResult, vector<string>> BooleanLPSolver::check(vector<Expression> cons
 		if (state().isBooleanVariable.at(index) || isConditionalConstraint(index))
 			resizeAndSet(booleanVariables, index, name);
 		else
-			resizeAndSet(lpState.variableNames, index, name);
+			lpSolver.setVariableName(index, name);
 
-	// TODO keep a cache as a member that is never reset.
-	// TODO We can also keep the split unconditionals across push/pop
-	// We only need to be careful to update the number of variables.
-
-	std::vector<std::pair<size_t, LPSolver>> lpSolvers;
-
-	// TODO We start afresh here. If we want this to reuse the existing results
-	// from previous invocations of the boolean solver, we still have to use
-	// a cache.
-	// The current optimization is only for CDCL.
-	lpSolvers.emplace_back(0, LPSolver{&m_lpCache});
-	if (
-		lpSolvers.back().second.setState(lpState) == LPResult::Infeasible ||
-		lpSolvers.back().second.check().first == LPResult::Infeasible
-	)
+	if (lpSolver.check().first == LPResult::Infeasible)
 	{
 		cout << "----->>>>> unsatisfiable" << endl;
 		return {CheckResult::UNSATISFIABLE, {}};
@@ -263,8 +263,7 @@ pair<CheckResult, vector<string>> BooleanLPSolver::check(vector<Expression> cons
 				continue;
 			// "reason" is already stored for those constraints.
 			Constraint const& constraint = state().conditionalConstraints.at(constraintIndex);
-			solAssert(constraint.reasons.size() == 1 && *constraint.reasons.begin() == constraintIndex);
-			lpSolvers.back().second.addConstraint(constraint);
+			lpSolvers.back().second.addConstraint(constraint, constraintIndex);
 		}
 		auto&& [result, modelOrReason] = lpSolvers.back().second.check();
 		// We can only really use the result "infeasible". Everything else should be "sat".
@@ -365,7 +364,7 @@ optional<Literal> BooleanLPSolver::parseLiteral(smtutil::Expression const& _expr
 		LinearExpression data = *left - *right;
 		data[0] *= -1;
 
-		return Literal{true, addConditionalConstraint(Constraint{move(data), _expr.name == "=", {}})};
+		return Literal{true, addConditionalConstraint(Constraint{move(data), _expr.name == "="})};
 	}
 	else if (_expr.name == ">=")
 		return parseLiteral(_expr.arguments.at(1) <= _expr.arguments.at(0));
@@ -390,7 +389,6 @@ Literal BooleanLPSolver::negate(Literal const& _lit)
 			Constraint le = c;
 			le.equality = false;
 			le.data[0] -= 1;
-			le.reasons.clear();
 			Literal leL{true, addConditionalConstraint(le)};
 
 			// X >= b + 1
@@ -399,7 +397,6 @@ Literal BooleanLPSolver::negate(Literal const& _lit)
 			ge.equality = false;
 			ge.data *= -1;
 			ge.data[0] -= 1;
-			ge.reasons.clear();
 			Literal geL{true, addConditionalConstraint(ge)};
 
 
@@ -419,7 +416,6 @@ Literal BooleanLPSolver::negate(Literal const& _lit)
 			Constraint negated = c;
 			negated.data *= -1;
 			negated.data[0] -= 1;
-			negated.reasons.clear();
 			return Literal{true, addConditionalConstraint(negated)};
 		}
 	}
@@ -561,8 +557,6 @@ size_t BooleanLPSolver::addConditionalConstraint(Constraint _constraint)
 	//  - integers
 	declareVariable(name, false);
 	size_t index = state().variables.at(name);
-	solAssert(_constraint.reasons.empty());
-	_constraint.reasons.emplace(index);
 	state().conditionalConstraints[index] = move(_constraint);
 	return index;
 }
diff --git a/libsolutil/BooleanLP.h b/libsolutil/BooleanLP.h
index 85e7efcde..282ad94ee 100644
--- a/libsolutil/BooleanLP.h
+++ b/libsolutil/BooleanLP.h
@@ -40,8 +40,14 @@ struct State
 	std::map<size_t, Constraint> conditionalConstraints;
 	std::vector<Clause> clauses;
 
+	struct Bounds
+	{
+		std::optional<rational> lower;
+		std::optional<rational> upper;
+	};
+
 	// Unconditional bounds on variables
-	std::map<size_t, SolvingState::Bounds> bounds;
+	std::map<size_t, Bounds> bounds;
 	// Unconditional constraints
 	std::vector<Constraint> fixedConstraints;
 };
@@ -122,9 +128,6 @@ private:
 
 	/// Stack of state, to allow for push()/pop().
 	std::vector<State> m_state{{State{}}};
-
-	std::unordered_map<SolvingState, LPResult> m_lpCache;
-
 };
 
 
diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp
index d8593667a..5173dffd0 100644
--- a/libsolutil/LP.cpp
+++ b/libsolutil/LP.cpp
@@ -63,17 +63,6 @@ inline std::vector<bool>& operator|=(std::vector<bool>& _x, std::vector<bool> co
 	return _x;
 }
 
-/**
- * Simplex tableau.
- */
-struct Tableau
-{
-	/// The factors of the objective function (first row of the tableau)
-	LinearExpression objective;
-	/// The tableau matrix (equational form).
-	std::vector<LinearExpression> data;
-};
-
 string toString(rational const& _x)
 {
 	if (_x == bigint(1) << 256)
@@ -95,414 +84,6 @@ string reasonToString(ReasonSet const& _reasons, size_t _minSize)
 	return result;
 }
 
-/*
-string toString(LinearExpression const& _expr)
-{
-	vector<string> items;
-	for (auto&& multiplier: _expr)
-		 if (multiplier != 0)
-			items.emplace_back(::toString(multiplier));
-		else
-			items.emplace_back("_");
-	for (string& item: items)
-		while (item.size() < 3)
-			item = " " + item;
-	return joinHumanReadable(items, " ");
-}
-
-string toString(Tableau const& _tableau)
-{
-	string s = toString(_tableau.objective) + "\n";
-	for (auto&& d: _tableau.data)
-		s += toString(d) + "\n";
-	return s;
-}
-
-
-/// Adds slack variables to remove non-equality costraints from a set of constraints
-/// and returns the data part of the tableau / constraints.
-/// The second return variable is true if the original input had any equality constraints.
-pair<vector<LinearExpression>, bool> toEquationalForm(vector<Constraint> _constraints)
-{
-	size_t varsNeeded = static_cast<size_t>(ranges::count_if(_constraints, [](Constraint const& _c) { return !_c.equality; }));
-	if (varsNeeded > 0)
-	{
-		size_t columns = _constraints.at(0).data.size();
-		size_t varsAdded = 0;
-		for (Constraint& constraint: _constraints)
-		{
-			solAssert(constraint.data.size() == columns, "");
-			constraint.data.resize(columns + varsNeeded);
-			if (!constraint.equality)
-			{
-				constraint.equality = true;
-				constraint.data[columns + varsAdded] = bigint(1);
-				varsAdded++;
-			}
-		}
-		solAssert(varsAdded == varsNeeded);
-	}
-
-	vector<LinearExpression> data;
-	for (Constraint& c: _constraints)
-		data.emplace_back(move(c.data));
-
-	return make_pair(move(data), varsNeeded < _constraints.size());
-}
-
-/// Finds the simplex pivot column: The column with the largest positive objective factor.
-/// If all objective factors are zero or negative, the optimum has been found and nullopt is returned.
-optional<size_t> findPivotColumn(Tableau const& _tableau)
-{
-	auto&& [maxColumn, maxValue] = ranges::max(
-		_tableau.objective.enumerateTail(),
-		{},
-		[](std::pair<size_t, rational> const& _x) { return _x.second; }
-	);
-
-	if (maxValue <= rational{0})
-		return nullopt; // found optimum
-	else
-		return maxColumn;
-}
-
-/// Finds the simplex pivot row, given the column:
-/// If there is no positive factor in the column, the problem is unbounded, nullopt is returned.
-/// Otherwise, returns the row i such that c[i] / x[i] is minimal and x[i] is positive, where
-/// c[i] is the constant factor (not the objective factor!) in row i.
-optional<size_t> findPivotRow(Tableau const& _tableau, size_t _pivotColumn)
-{
-	auto positiveColumnEntries =
-		ranges::views::iota(size_t(0), _tableau.data.size()) |
-		ranges::views::transform([&](size_t i) {
-			return make_pair(i, _tableau.data[i][_pivotColumn]);
-		}) |
-		ranges::views::filter([](pair<size_t, rational> const& _entry) {
-			return _entry.second.numerator() > 0;
-		});
-	if (positiveColumnEntries.empty())
-		return nullopt; // unbounded
-
-	return ranges::min(
-		positiveColumnEntries,
-		{},
-		[&](std::pair<size_t, rational> const& _entry) {
-			return _tableau.data[_entry.first][0] / _entry.second;
-		}
-	).first;
-}
-
-/// Performs equivalence transform on @a _tableau, so that
-/// the column @a _pivotColumn is all zeros (including the objective row) except for @a _pivotRow,
-/// where it is 1.
-void performPivot(Tableau& _tableau, size_t _pivotRow, size_t _pivotColumn)
-{
-	rational pivot = _tableau.data[_pivotRow][_pivotColumn];
-	solAssert(pivot != 0, "");
-	if (pivot != 1)
-		_tableau.data[_pivotRow] /= pivot;
-	solAssert(_tableau.data[_pivotRow][_pivotColumn] == rational(1), "");
-
-	LinearExpression const& _pivotRowData = _tableau.data[_pivotRow];
-	auto subtractMultipleOfPivotRow = [&](LinearExpression& _row) {
-		if (_row[_pivotColumn] == rational{1})
-			_row -= _pivotRowData;
-		else if (_row[_pivotColumn] == rational{-1})
-			_row += _pivotRowData;
-		else if (_row[_pivotColumn])
-			_row -= _row[_pivotColumn] * _pivotRowData;
-	};
-
-	subtractMultipleOfPivotRow(_tableau.objective);
-	for (size_t i = 0; i < _tableau.data.size(); ++i)
-		if (i != _pivotRow)
-			subtractMultipleOfPivotRow(_tableau.data[i]);
-}
-
-/// Transforms the tableau such that the last vectors are basis vectors
-/// and their objective coefficients are zero.
-/// Makes various assumptions and should only be used after adding
-/// a certain number of slack variables.
-void selectLastVectorsAsBasis(Tableau& _tableau)
-{
-	// We might skip the operation for a column if it is already the correct
-	// unit vector and its objective coefficient is zero.
-	size_t columns = _tableau.objective.size();
-	size_t rows = _tableau.data.size();
-	for (size_t i = 0; i < rows; ++i)
-		performPivot(_tableau, i, columns - rows + i);
-}
-
-/// If column @a _column inside tableau is a basis vector
-/// (i.e. one entry is 1, the others are 0), returns the index
-/// of the 1, otherwise nullopt.
-optional<size_t> basisIndex(Tableau const& _tableau, size_t _column)
-{
-	optional<size_t> row;
-	for (size_t i = 0; i < _tableau.data.size(); ++i)
-		if (_tableau.data[i][_column] == bigint(1))
-		{
-			if (row)
-				return std::nullopt;
-			else
-				row = i;
-		}
-		else if (_tableau.data[i][_column] != 0)
-			return std::nullopt;
-	return row;
-}
-
-/// @returns a solution vector, assuming one exists.
-/// The solution vector minimizes the objective function if the tableau
-/// is the result of the simplex algorithm.
-vector<rational> solutionVector(Tableau const& _tableau)
-{
-	vector<rational> result;
-	vector<bool> rowsSeen(_tableau.data.size(), false);
-	for (size_t j = 1; j < _tableau.objective.size(); j++)
-	{
-		optional<size_t> row = basisIndex(_tableau, j);
-		if (row && rowsSeen[*row])
-			row = nullopt;
-		result.emplace_back(row ? _tableau.data[*row][0] : rational{});
-		if (row)
-			rowsSeen[*row] = true;
-	}
-	return result;
-}
-
-
-/// Solve the LP A x = b s.t. min c^T x
-/// Here, c is _tableau.objective and the first column of _tableau.data
-/// encodes b and the other columns encode A
-/// Assumes the tableau has a trivial basic feasible solution.
-/// Tries for a number of iterations and then gives up.
-pair<LPResult, Tableau> simplexEq(Tableau _tableau)
-{
-	size_t const iterations = min<size_t>(60, 50 + _tableau.objective.size() * 2);
-	for (size_t step = 0; step <= iterations; ++step)
-	{
-		optional<size_t> pivotColumn = findPivotColumn(_tableau);
-		if (!pivotColumn)
-			return make_pair(LPResult::Feasible, move(_tableau));
-		optional<size_t> pivotRow = findPivotRow(_tableau, *pivotColumn);
-		if (!pivotRow)
-			return make_pair(LPResult::Unbounded, move(_tableau));
-		performPivot(_tableau, *pivotRow, *pivotColumn);
-	}
-	return make_pair(LPResult::Unknown, Tableau{});
-}
-
-/// We add slack variables to find a basic feasible solution.
-/// In particular, there is a slack variable for each row
-/// which is weighted negatively. Setting the new slack
-/// variables to one and all other variables to zero yields
-/// a basic feasible solution.
-/// If the optimal solution has all slack variables set to zero,
-/// this is a basic feasible solution. Otherwise, the original
-/// problem is infeasible.
-/// This function returns the modified tableau with the original
-/// objective function and the slack variables removed.
-pair<LPResult, Tableau> simplexPhaseI(Tableau _tableau)
-{
-	LinearExpression originalObjective = _tableau.objective;
-
-	size_t rows = _tableau.data.size();
-	size_t columns = _tableau.objective.size();
-	for (size_t i = 0; i < rows; ++i)
-	{
-		if (_tableau.data[i][0] < 0)
-			_tableau.data[i] *= -1;
-		_tableau.data[i].resize(columns + rows);
-		_tableau.data[i][columns + i] = 1;
-	}
-	_tableau.objective = {};
-	_tableau.objective.resize(columns);
-	_tableau.objective.resize(columns + rows, rational{-1});
-
-	// This sets the objective factors of the slack variables
-	// to zero (and thus selects a basic feasible solution).
-	selectLastVectorsAsBasis(_tableau);
-
-	LPResult result;
-	tie(result, _tableau) = simplexEq(move(_tableau));
-	if (result == LPResult::Unknown)
-		return make_pair(LPResult::Unknown, Tableau{});
-	solAssert(result != LPResult::Infeasible, "");
-
-	vector<rational> optimum = solutionVector(_tableau);
-
-	// If the solution needs a nonzero factor for a slack variable,
-	// the original system is infeasible.
-	for (size_t i = columns - 1; i < optimum.size(); ++i)
-		if (optimum[i] != 0)
-			return make_pair(LPResult::Infeasible, Tableau{});
-
-	// Restore original objective and remove slack variables.
-	_tableau.objective = move(originalObjective);
-	for (auto& row: _tableau.data)
-		row.resize(columns);
-
-	return make_pair(LPResult::Feasible, move(_tableau));
-}
-
-/// Returns true if the all-zero solution is not a solution for the tableau.
-bool needsPhaseI(Tableau const& _tableau)
-{
-	for (auto const& row: _tableau.data)
-		if (row[0] < 0)
-			return true;
-	return false;
-}
-
-
-/// Solve the LP Ax <= b s.t. min c^Tx
-pair<LPResult, vector<rational>> simplex(vector<Constraint> _constraints, LinearExpression _objectives)
-{
-	Tableau tableau;
-	tableau.objective = move(_objectives);
-	bool hasEquations = false;
-	tie(tableau.data, hasEquations) = toEquationalForm(_constraints);
-	tableau.objective.resize(tableau.data.at(0).size());
-	//cout << "Running simplex on " << tableau.objective.size() << " x " << tableau.data.size() << endl;
-
-	if (hasEquations || needsPhaseI(tableau))
-	{
-		LPResult result;
-		tie(result, tableau) = simplexPhaseI(move(tableau));
-		if (result == LPResult::Infeasible || result == LPResult::Unknown)
-			return make_pair(result, vector<rational>{});
-		solAssert(result == LPResult::Feasible, "");
-	}
-	// We know that the system is satisfiable and we know a solution,
-	// but it might not be optimal.
-	LPResult result;
-	tie(result, tableau) = simplexEq(move(tableau));
-	solAssert(result != LPResult::Infeasible, "");
-	return make_pair(result, solutionVector(tableau));
-}
-
-
-/// Turns all bounds into constraints.
-/// @returns false if the bounds make the state infeasible.
-optional<ReasonSet> boundsToConstraints(SolvingState& _state)
-{
-	size_t columns = _state.variableNames.size();
-
-	// Bound zero should not exist because the variable zero does not exist.
-	for (auto const& [varIndex, bounds]: _state.bounds | ranges::views::enumerate | ranges::views::tail)
-	{
-		if (bounds.lower && bounds.upper)
-		{
-			if (*bounds.lower > *bounds.upper)
-				return bounds.lowerReasons + bounds.upperReasons;
-			if (*bounds.lower == *bounds.upper)
-			{
-				LinearExpression c;
-				c.resize(columns);
-				c[0] = *bounds.lower;
-				c[varIndex] = bigint(1);
-				_state.constraints.emplace_back(Constraint{move(c), true, bounds.lowerReasons + bounds.upperReasons});
-				continue;
-			}
-		}
-		if (bounds.lower && *bounds.lower > 0)
-		{
-			LinearExpression c;
-			c.resize(columns);
-			c[0] = -*bounds.lower;
-			c[varIndex] = bigint(-1);
-			_state.constraints.emplace_back(Constraint{move(c), false, move(bounds.lowerReasons)});
-		}
-		if (bounds.upper)
-		{
-			LinearExpression c;
-			c.resize(columns);
-			c[0] = *bounds.upper;
-			c[varIndex] = bigint(1);
-			_state.constraints.emplace_back(Constraint{move(c), false, move(bounds.upperReasons)});
-		}
-	}
-	_state.bounds.clear();
-	return nullopt;
-}
-*/
-
-/// Removes incides set to true from a vector-like data structure.
-template <class T>
-void eraseIndices(T& _data, vector<bool> const& _indicesToRemove)
-{
-	T result;
-	for (size_t i = 0; i < _data.size(); i++)
-		if (!_indicesToRemove[i])
-			result.push_back(move(_data[i]));
-	_data = move(result);
-}
-
-
-void removeColumns(SolvingState& _state, vector<bool> const& _columnsToRemove)
-{
-	eraseIndices(_state.bounds, _columnsToRemove);
-	for (Constraint& constraint: _state.constraints)
-		eraseIndices(constraint.data, _columnsToRemove);
-	eraseIndices(_state.variableNames, _columnsToRemove);
-}
-
-auto nonZeroEntriesInColumn(SolvingState const& _state, size_t _column)
-{
-	return
-		_state.constraints |
-		ranges::views::enumerate |
-		ranges::views::filter([=](auto const& _entry) { return _entry.second.data[_column]; }) |
-		ranges::views::transform([](auto const& _entry) { return _entry.first; });
-}
-
-/// @returns vectors of column- and row-indices that are connected to the given column,
-/// in the sense of variables occurring in a constraint and constraints for variables.
-pair<vector<bool>, vector<bool>> connectedComponent(SolvingState const& _state, size_t _column)
-{
-	solAssert(_state.variableNames.size() >= 2, "");
-
-	vector<bool> includedColumns(_state.variableNames.size(), false);
-	vector<bool> seenColumns(_state.variableNames.size(), false);
-	vector<bool> includedRows(_state.constraints.size(), false);
-	stack<size_t> columnsToProcess;
-	columnsToProcess.push(_column);
-	while (!columnsToProcess.empty())
-	{
-		size_t column = columnsToProcess.top();
-		columnsToProcess.pop();
-		if (includedColumns[column])
-			continue;
-		includedColumns[column] = true;
-
-		for (size_t row: nonZeroEntriesInColumn(_state, column))
-		{
-			if (includedRows[row])
-				continue;
-			includedRows[row] = true;
-			for (auto const& [index, entry]: _state.constraints[row].data.enumerateTail())
-				if (entry && !seenColumns[index])
-				{
-					seenColumns[index] = true;
-					columnsToProcess.push(index);
-				}
-		}
-	}
-	return make_pair(move(includedColumns), move(includedRows));
-}
-
-void normalizeRowLengths(SolvingState& _state)
-{
-	size_t vars = max(_state.variableNames.size(), _state.bounds.size());
-	for (Constraint const& c: _state.constraints)
-		vars = max(vars, c.data.size());
-	_state.variableNames.resize(vars);
-	_state.bounds.resize(vars);
-	for (Constraint& c: _state.constraints)
-		c.data.resize(vars);
-}
-
 }
 
 
@@ -556,8 +137,6 @@ set<size_t> SolvingState::reasons() const
 	set<size_t> ret;
 	for (Bounds const& b: bounds)
 		ret += b.lowerReasons + b.upperReasons;
-	for (Constraint const& c: constraints)
-		ret += c.reasons;
 	return ret;
 }
 
@@ -580,7 +159,6 @@ string SolvingState::toString() const
 				line.emplace_back(mult + variableNames.at(index));
 			}
 		result +=
-			reasonToString(constraint.reasons, reasonLength) +
 			joinHumanReadable(line, " + ") +
 			(constraint.equality ? "  = " : " <= ") +
 			::toString(constraint.data.front()) +
@@ -603,181 +181,429 @@ string SolvingState::toString() const
 	return result;
 }
 
-SimplexWithBounds::SimplexWithBounds(SolvingState _state):
-	m_state(move(_state))
+
+void LPSolver::addConstraint(Constraint const& _constraint, optional<size_t> _reason)
 {
-	size_t varsNeeded = m_state.constraints.size();
-	size_t row = 0;
-	// TODO the new method does not need zero as lower-bound on all variables,
-	// but we currently assume that in the BooleanLP at some places.
-	for (auto& bounds: m_state.bounds)
-		if (!bounds.lower)
-			bounds.lower = rational{0};
+	// TODO at this point, we could also determine if it is a fixed variable.
+	// (maybe even taking the bounds on existing variables into account)
+	// If we do this, we have to take reasons into account properly!
+	set<size_t> touchedProblems;
+	for (auto const& [index, entry]: _constraint.data.enumerateTail())
+		if (entry)
+			if (m_subProblemsPerVariable.count(index))
+				touchedProblems.emplace(m_subProblemsPerVariable.at(index));
 
-	for (Constraint& c: m_state.constraints)
+	if (touchedProblems.empty())
 	{
-		// TODO duplicated in addConstraint
-		c.data.resize(m_state.variableNames.size() + varsNeeded);
-		size_t newVarIndex = m_state.variableNames.size();
-		m_basicVariables[newVarIndex] = row++;
-		solAssert(m_state.variableNames.size() == m_state.bounds.size());
-		// TODO name needed unique?
-		m_state.variableNames.emplace_back("_s" + to_string(newVarIndex));
-		m_state.bounds.emplace_back();
-		if (c.equality)
-			m_state.bounds.back().lower = c.data[0];
-		m_state.bounds.back().upper = c.data[0];
-		c.equality = true;
-		solAssert(c.data.size() > newVarIndex);
-		c.data[newVarIndex] = -1;
-		c.data[0] = 0;
+		//cout << "Creating new sub problem." << endl;
+		// TODO we could find an empty spot for the pointer.
+		m_subProblems.emplace_back(make_shared<SubProblem>());
+		solAssert(!m_subProblems.back()->sealed);
+		touchedProblems.emplace(m_subProblems.size() - 1);
 	}
-
-	m_assignments.resize(m_state.variableNames.size());
-	solAssert(m_assignments.size() == m_state.bounds.size());
+	for (size_t problemToErase: touchedProblems | ranges::views::tail | ranges::views::reverse)
+		combineSubProblems(*touchedProblems.begin(), problemToErase);
+	addConstraintToSubProblem(*touchedProblems.begin(), _constraint, move(_reason));
+	//cout << "Added constraint:\n" << toString() << endl;
 }
 
-LPResult SimplexWithBounds::check()
+void LPSolver::setVariableName(size_t _variable, string _name)
 {
-	solAssert(m_assignments.size() == m_state.bounds.size());
-	cout << "===========================" << endl;
-	cout << "===========================" << endl;
-	cout << "simpl with bounds on\n" << toString() << endl;
+	// TODO it might be constly to do this before we know hich variables relate to which
 
-	// We start with an all-zero assignment and then gradually add
-	// the bounds we already have.
+	SubProblem& p = unsealForVariable(_variable);
+	p.variables[p.varMapping.at(_variable)].name = move(_name);
+}
 
-	//cout << "Adjusting bounds on non-basic variables\n";
-	// Adjust the assignments so we satisfy the bounds of the non-basic variables.
-	for (auto const& [i, bounds]: m_state.bounds | ranges::views::enumerate)
+void LPSolver::addLowerBound(size_t _variable, rational _bound)
+{
+	SubProblem& p = unsealForVariable(_variable);
+	Variable& var = p.variables[p.varMapping.at(_variable)];
+	if (!var.bounds.lower || *var.bounds.lower < _bound)
+		var.bounds.lower = move(_bound);
+}
+
+void LPSolver::addUpperBound(size_t _variable, rational _bound)
+{
+	SubProblem& p = unsealForVariable(_variable);
+	Variable& var = p.variables[p.varMapping.at(_variable)];
+	if (!var.bounds.upper || *var.bounds.upper > _bound)
+		var.bounds.upper = move(_bound);
+}
+
+pair<LPResult, variant<Model, ReasonSet>> LPSolver::check()
+{
+	for (auto&& [index, problem]: m_subProblems | ranges::views::enumerate)
+		if (problem)
+			problem->sealed = true;
+
+	for (auto&& [index, problem]: m_subProblems | ranges::views::enumerate)
 	{
-		if (m_basicVariables.count(i) || (!bounds.lower && !bounds.upper))
+		if (!problem)
 			continue;
-		if (bounds.lower && bounds.upper)
-			solAssert(*bounds.lower <= *bounds.upper);
-		if (bounds.lower && m_assignments[i] < *bounds.lower)
-			update(i, *bounds.lower);
-		else if (bounds.upper && m_assignments[i] > *bounds.upper)
-			update(i, *bounds.upper);
+		if (!problem->result)
+			problem->result = problem->check();
+
+		if (*problem->result == LPResult::Infeasible)
+			return {LPResult::Infeasible, problem->reasons};
+	}
+	//cout << "Feasible:\n" << toString() << endl;
+	return {LPResult::Feasible, model()};
+}
+
+string LPSolver::toString() const
+{
+	string result = "LP Solver state:\n";
+	for (auto const& problem: m_subProblems)
+		if (problem)
+			result += problem->toString();
+	return result;
+}
+
+LPSolver::SubProblem& LPSolver::unseal(size_t _problemIndex)
+{
+	shared_ptr<SubProblem>& problem = m_subProblems[_problemIndex];
+	solAssert(problem);
+	if (problem->sealed)
+		problem = make_shared<SubProblem>(*problem);
+	problem->sealed = false;
+	problem->result = nullopt;
+	return *problem;
+}
+
+LPSolver::SubProblem& LPSolver::unsealForVariable(size_t _outerIndex)
+{
+	if (!m_subProblemsPerVariable.count(_outerIndex))
+	{
+		m_subProblems.emplace_back(make_shared<SubProblem>());
+		addOuterVariableToSubProblem(m_subProblems.size() - 1, _outerIndex);
+	}
+	return unseal(m_subProblemsPerVariable.at(_outerIndex));
+}
+
+void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom)
+{
+	//cout << "Combining\n" << m_subProblems.at(_combineFrom)->toString();
+	//cout << "\ninto\n" << m_subProblems.at(_combineInto)->toString();
+	SubProblem& combineInto = unseal(_combineInto);
+	SubProblem const& combineFrom = *m_subProblems[_combineFrom];
+
+	size_t varShift = combineInto.variables.size();
+	size_t rowShift = combineInto.factors.size();
+	size_t newRowLength = combineInto.variables.size() + combineFrom.variables.size();
+	for (LinearExpression& row: combineInto.factors)
+		row.resize(newRowLength);
+	for (LinearExpression const& row: combineFrom.factors)
+	{
+		LinearExpression shiftedRow;
+		shiftedRow.resize(newRowLength);
+		for (auto&& [index, f]: row.enumerate())
+			shiftedRow[varShift + index] = f;
+		combineInto.factors.emplace_back(move(shiftedRow));
+	}
+	combineInto.variables += combineFrom.variables;
+	for (auto&& [index, row]: combineFrom.basicVariables)
+		combineInto.basicVariables.emplace(index + varShift, row + rowShift);
+	for (auto&& [outerIndex, innerIndex]: combineFrom.varMapping)
+		combineInto.varMapping.emplace(outerIndex, innerIndex + varShift);
+	combineInto.reasons += combineFrom.reasons;
+
+	for (auto& item: m_subProblemsPerVariable)
+		if (item.second == _combineFrom)
+			item.second = _combineInto;
+
+	m_subProblems[_combineFrom].reset();
+	//cout << "result: \n" << m_subProblems.at(_combineInto)->toString();
+	//cout << "------------------------------\n";
+}
+
+// TODO move this function into the problem struct and make it erturn set of vaiables added
+
+void LPSolver::addConstraintToSubProblem(
+	size_t _subProblem,
+	Constraint const& _constraint,
+	std::optional<size_t> _reason
+)
+{
+	// TODO opt:
+	// Add "fixed variables" at general state (above sub problems)
+	// replace all fixed variables in constraint by their values
+	// If constaint is direct constraint on variable, just add it to its bounds
+	// If constraint results in variable being fixed,
+	// then push that to the general state as 'fixed variables'
+	// Remove the variable from the subproblem (so we can more efficiently split)
+	// If we remove the var, it is a bit more tricky because we have to store the reason
+	// together with the var.
+
+	SubProblem& problem = unseal(_subProblem);
+	if (_reason)
+		problem.reasons.insert(*_reason);
+
+	size_t numVariables = 0;
+	size_t latestVariableIndex = size_t(-1);
+	// Make all variables available and check if it is a simple bound on a variable.
+	for (auto const& [index, entry]: _constraint.data.enumerateTail())
+		if (entry)
+		{
+			latestVariableIndex = index;
+			numVariables++;
+			if (!problem.varMapping.count(index))
+				addOuterVariableToSubProblem(_subProblem, index);
+		}
+	if (numVariables == 1)
+	{
+		// Add this as direct bound.
+		// TODO we could avoid some of the steps by introducing an "addUpperBound"
+		// function on the subproblem.
+		rational factor = _constraint.data[latestVariableIndex];
+		rational bound = _constraint.data.front() / factor;
+		if (factor > 0 || _constraint.equality)
+			addUpperBound(latestVariableIndex, bound);
+		if (factor < 0 || _constraint.equality)
+			addLowerBound(latestVariableIndex, bound);
+		return;
+	}
+
+	// Introduce the slack variable.
+	size_t slackIndex = addNewVariableToSubProblem(_subProblem);
+	// Name is only needed for printing
+	//problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++);
+	problem.basicVariables[slackIndex] = problem.factors.size();
+	if (_constraint.equality)
+		problem.variables[slackIndex].bounds.lower = _constraint.data[0];
+	problem.variables[slackIndex].bounds.upper = _constraint.data[0];
+
+
+	// Compress the constraint, i.e. turn outer variable indices into
+	// inner variable indices.
+	rational valueForSlack;
+	LinearExpression compressedConstraint;
+	LinearExpression basicVarNullifier;
+	compressedConstraint.resize(problem.variables.size());
+	for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail())
+		if (entry)
+		{
+			size_t innerIndex = problem.varMapping.at(outerIndex);
+			if (problem.basicVariables.count(innerIndex))
+			{
+				// We cannot add basic variables directly, so replace them by their row.
+				basicVarNullifier += entry * problem.factors.at(problem.basicVariables[innerIndex]);
+				basicVarNullifier[innerIndex] = {};
+			}
+			else
+				compressedConstraint[innerIndex] = entry;
+			valueForSlack += problem.variables[innerIndex].value * entry;
+		}
+	compressedConstraint += move(basicVarNullifier);
+
+	compressedConstraint[slackIndex] = -1;
+	problem.factors.emplace_back(move(compressedConstraint));
+	problem.basicVariables[slackIndex] = problem.factors.size() - 1;
+	problem.variables[slackIndex].value = valueForSlack;
+}
+
+void LPSolver::addOuterVariableToSubProblem(size_t _subProblem, size_t _outerIndex)
+{
+	size_t index = addNewVariableToSubProblem(_subProblem);
+	unseal(_subProblem).varMapping.emplace(_outerIndex, index);
+	m_subProblemsPerVariable[_outerIndex] = _subProblem;
+}
+
+size_t LPSolver::addNewVariableToSubProblem(size_t _subProblem)
+{
+	SubProblem& problem = unseal(_subProblem);
+	size_t index = problem.variables.size();
+	for (LinearExpression& c: problem.factors)
+		c.resize(index + 1);
+	problem.variables.emplace_back();
+	return index;
+}
+
+map<string, rational> LPSolver::model() const
+{
+	map<string, rational> result;
+	/*
+	 * // This is actually unused
+	for (auto const& problem: m_subProblems)
+		if (problem)
+			for (auto&& [outerIndex, innerIndex]: problem->varMapping)
+				result[problem->variables[innerIndex].name] = problem->variables[innerIndex].value;
+	cout << "MOdel size: " << to_string(result.size()) << endl;
+	*/
+	return result;
+}
+
+LPResult LPSolver::SubProblem::check()
+{
+
+	// TODO one third of the computing time (inclusive) in this function
+	// is spent on "operator<" - maybe we can cache "is in bounds" for variables
+	// and invalidate that in the update procedures.
+
+//	cout << "checking..." << endl;
+//	cout << toString() << endl;
+//	cout << "----------------------------" << endl;
+//	cout << "fixing non-basic..." << endl;
+	// Adjust the assignments so we satisfy the bounds of the non-basic variables.
+	for (auto const& [i, var]: variables | ranges::views::enumerate)
+	{
+		if (var.bounds.lower && var.bounds.upper && *var.bounds.lower > *var.bounds.upper)
+			return LPResult::Infeasible;
+		if (basicVariables.count(i) || (!var.bounds.lower && !var.bounds.upper))
+			continue;
+		if (var.bounds.lower && var.value < *var.bounds.lower)
+			update(i, *var.bounds.lower);
+		else if (var.bounds.upper && var.value > *var.bounds.upper)
+			update(i, *var.bounds.upper);
 	}
-	//cout << "Bounds on non-basic vaiables set.\n" << toString() << endl;
 
 	// Now try to make the basic variables happy, pivoting if necessary.
 
+//	cout << "fixed non-basic." << endl;
+//	cout << toString() << endl;
+//	cout << "----------------------------" << endl;
+
 	// TODO bound number of iterations
 	while (auto bvi = firstConflictingBasicVariable())
 	{
-		//cout << "Basic variable " << m_state.variableNames[*bvi] << " is conflicting." << endl;
-		if (m_state.bounds[*bvi].lower && m_assignments[*bvi] < *m_state.bounds[*bvi].lower)
+		Variable const& basicVar = variables[*bvi];
+//		cout << toString() << endl;
+//		cout << "Fixing basic " << basicVar.name << endl;
+//		cout << "----------------------------" << endl;
+		if (basicVar.bounds.lower && basicVar.bounds.upper)
+			solAssert(*basicVar.bounds.lower <= *basicVar.bounds.upper);
+		if (basicVar.bounds.lower && basicVar.value < *basicVar.bounds.lower)
 		{
 			if (auto replacementVar = firstReplacementVar(*bvi, true))
-				pivotAndUpdate(*bvi, *m_state.bounds.at(*bvi).lower, *replacementVar);
-			else
 			{
-				//cout << "------->>>>>>>>>>>>>> infeasible\n";
-				return LPResult::Infeasible;
+//				cout << "Replacing by " << variables[*replacementVar].name << endl;
+
+				pivotAndUpdate(*bvi, *basicVar.bounds.lower, *replacementVar);
 			}
+			else
+				return LPResult::Infeasible;
 		}
-		else if (m_state.bounds[*bvi].upper && m_assignments[*bvi] > *m_state.bounds[*bvi].upper)
+		else if (basicVar.bounds.upper && basicVar.value > *basicVar.bounds.upper)
 		{
 			if (auto replacementVar = firstReplacementVar(*bvi, false))
-				pivotAndUpdate(*bvi, *m_state.bounds.at(*bvi).upper, *replacementVar);
-			else
 			{
-				//cout << "------->>>>>>>>>>>>>> infeasible\n";
-				return LPResult::Infeasible;
+//				cout << "Replacing by " << variables[*replacementVar].name << endl;
+				pivotAndUpdate(*bvi, *basicVar.bounds.upper, *replacementVar);
 			}
+			else
+				return LPResult::Infeasible;
 		}
+//		cout << "Fixed basic " << basicVar.name << endl;
+//		cout << toString() << endl;
 	}
 
-	//cout << "SAT ------->>>>>>>>>>>>>>\n";
-
 	return LPResult::Feasible;
+}
+
+string LPSolver::SubProblem::toString() const
+{
+	string resultString;
+	for (auto&& [i, v]: variables | ranges::views::enumerate)
+	{
+		if (v.bounds.lower)
+			resultString += ::toString(*v.bounds.lower) + " <= ";
+		else
+			resultString += "       ";
+		resultString += v.name;
+		if (v.bounds.upper)
+			resultString += " <= " + ::toString(*v.bounds.upper);
+		else
+			resultString += "       ";
+		resultString += "   := " + ::toString(v.value) + "\n";
+	}
+	for (auto&& [rowIndex, row]: factors | ranges::views::enumerate)
+	{
+		string basicVarPrefix;
+		string rowString;
+		for (auto&& [i, f]: row.enumerate())
+		{
+			if (basicVariables.count(i) && basicVariables.at(i) == rowIndex)
+			{
+				solAssert(f == -1);
+				solAssert(basicVarPrefix.empty());
+				basicVarPrefix = variables[i].name + " = ";
+			}
+			else if (f != 0)
+			{
+				string joiner = f < 0 ? " - " : f > 0 && !rowString.empty() ? " + " : " ";
+				string factor = f == 1 || f == -1 ? "" : ::toString(abs(f)) + " ";
+				string var = variables[i].name;
+				rowString += joiner + factor + var;
+			}
+		}
+		resultString += basicVarPrefix + rowString + "\n";
+	}
+	if (result)
+	{
+		if (*result == LPResult::Feasible)
+			resultString += "result: feasible\n";
+		else
+			resultString += "result: infeasible\n";
+	}
+	else
+		resultString += "result: unknown\n";
+
+
+	return resultString + "----\n";
+}
+
+void LPSolver::SubProblem::update(size_t _varIndex, rational const& _value)
+{
+	rational delta = _value - variables[_varIndex].value;
+	variables[_varIndex].value = _value;
+	for (size_t j = 0; j < variables.size(); j++)
+		if (basicVariables.count(j) && factors[basicVariables.at(j)][_varIndex])
+			variables[j].value += delta * factors[basicVariables.at(j)][_varIndex];
 
 }
 
-size_t SimplexWithBounds::addVariable(string _name)
+optional<size_t> LPSolver::SubProblem::firstConflictingBasicVariable() const
 {
-	m_state.variableNames.emplace_back(move(_name));
-	m_state.bounds.emplace_back();
-	solAssert(m_state.variableNames.size() == m_state.bounds.size());
-	// TODO relax in the future.
-	m_state.bounds.back().lower = 0;
-	m_assignments.emplace_back(0);
-	size_t varsBefore = m_state.variableNames.size();
-	normalizeRowLengths(m_state);
-	solAssert(m_state.variableNames.size() == varsBefore);
-	m_assignments.resize(m_state.variableNames.size());
-	solAssert(m_assignments.size() == m_state.bounds.size());
-	return m_state.variableNames.size() - 1;
-}
-
-void SimplexWithBounds::addConstraint(Constraint _constraint)
-{
-	m_state.constraints.emplace_back(move(_constraint));
-	Constraint& c = m_state.constraints.back();
-	// TODO name needed unique?
-	size_t slack = addVariable("_s" + to_string(m_state.variableNames.size()));
-	c.data.resize(m_state.variableNames.size());
-	m_basicVariables[slack] = m_state.constraints.size() - 1;
-	if (c.equality)
-		m_state.bounds.at(slack).lower = c.data[0];
-	m_state.bounds.at(slack).upper = c.data[0];
-	c.equality = true;
-	solAssert(c.data.size() > slack);
-	c.data[slack] = -1;
-	c.data[0] = 0;
-	solAssert(m_assignments.size() == m_state.bounds.size());
-}
-
-void SimplexWithBounds::update(size_t _var, rational const& _value)
-{
-	rational delta = _value - m_assignments[_var];
-	m_assignments[_var] = _value;
-	for (size_t j = 0; j < m_assignments.size(); j++)
-		if (m_basicVariables.count(j) && m_state.constraints[m_basicVariables.at(j)].data[_var])
-			m_assignments[j] += m_state.constraints[m_basicVariables.at(j)].data[_var] * delta;
-
-}
-
-optional<size_t> SimplexWithBounds::firstConflictingBasicVariable() const
-{
-	for (auto const& [i, bounds]: m_state.bounds | ranges::views::enumerate)
-		if (m_basicVariables.count(i) && (
-			(bounds.lower && m_assignments[i] < *bounds.lower) ||
-			(bounds.upper && m_assignments[i] > *bounds.upper)
-		))
-			return i;
+	for (auto const& varItem: basicVariables)
+	{
+		Variable const& variable = variables[varItem.first];
+		if (
+			(variable.bounds.lower && variable.value < *variable.bounds.lower) ||
+			(variable.bounds.upper && variable.value > *variable.bounds.upper)
+		)
+			return varItem.first;
+	}
 	return nullopt;
 }
 
-optional<size_t> SimplexWithBounds::firstReplacementVar(
+optional<size_t> LPSolver::SubProblem::firstReplacementVar(
 	size_t _basicVarToReplace,
 	bool _increasing
 ) const
 {
-	LinearExpression const& basicVarEquation = m_state.constraints[m_basicVariables.at(_basicVarToReplace)].data;
-	for (auto const& [i, bounds]: m_state.bounds | ranges::views::enumerate)
+	LinearExpression const& basicVarEquation = factors[basicVariables.at(_basicVarToReplace)];
+	for (auto const& [i, var]: variables | ranges::views::enumerate)
 	{
-		if (m_basicVariables.count(i) || !basicVarEquation[i])
+		if (basicVariables.count(i) || !basicVarEquation[i])
 			continue;
 		bool positive = basicVarEquation[i] > 0;
 		if (!_increasing)
 			positive = !positive;
-		if (positive && (!m_state.bounds[i].upper || m_assignments[i] < m_state.bounds[i].upper))
+		Variable const& candidate = variables[i];
+		if (positive && (!candidate.bounds.upper || candidate.value < *candidate.bounds.upper))
 			return i;
-		if (!positive && (!m_state.bounds[i].lower || m_assignments[i] > m_state.bounds[i].lower))
+		if (!positive && (!candidate.bounds.lower || candidate.value > *candidate.bounds.lower))
 			return i;
 	}
 	return nullopt;
 }
 
-void SimplexWithBounds::pivot(size_t _old, size_t _new)
+void LPSolver::SubProblem::pivot(size_t _old, size_t _new)
 {
 	// Transform pivotRow such that the coefficient for _new is -1
 	// Then use that to set all other coefficients for _new to zero.
-	size_t pivotRow = m_basicVariables[_old];
-	LinearExpression& pivotRowData = m_state.constraints[m_basicVariables[_old]].data;
+	size_t pivotRow = basicVariables[_old];
+	LinearExpression& pivotRowData = factors[pivotRow];
 
 	rational pivot = pivotRowData[_new];
 	solAssert(pivot != 0, "");
@@ -796,572 +622,28 @@ void SimplexWithBounds::pivot(size_t _old, size_t _new)
 			_row += _row[_new] * pivotRowData;
 	};
 
-	for (size_t i = 0; i < m_state.constraints.size(); ++i)
+	for (size_t i = 0; i < factors.size(); ++i)
 		if (i != pivotRow)
-			subtractMultipleOfPivotRow(m_state.constraints[i].data);
+			subtractMultipleOfPivotRow(factors[i]);
 
-	m_basicVariables.erase(_old);
-	m_basicVariables[_new] = pivotRow;
+	basicVariables.erase(_old);
+	basicVariables[_new] = pivotRow;
 }
 
-void SimplexWithBounds::pivotAndUpdate(
+void LPSolver::SubProblem::pivotAndUpdate(
 	size_t _oldBasicVar,
 	rational const& _newValue,
 	size_t _newBasicVar
 )
 {
-	rational theta = (_newValue - m_assignments[_oldBasicVar]) / m_state.constraints[m_basicVariables[_oldBasicVar]].data[_newBasicVar];
+	rational theta = (_newValue - variables[_oldBasicVar].value) / factors[basicVariables[_oldBasicVar]][_newBasicVar];
 
-	m_assignments[_oldBasicVar] = _newValue;
-	m_assignments[_newBasicVar] += theta;
+	variables[_oldBasicVar].value = _newValue;
+	variables[_newBasicVar].value += theta;
 
-	for (auto const& [i, row]: m_basicVariables)
-		if (i != _oldBasicVar && m_state.constraints[row].data[_newBasicVar])
-			m_assignments[i] += m_state.constraints[row].data[_newBasicVar] * theta;
+	for (auto const& [i, row]: basicVariables)
+		if (i != _oldBasicVar && factors[row][_newBasicVar])
+			variables[i].value += factors[row][_newBasicVar] * theta;
 
 	pivot(_oldBasicVar, _newBasicVar);
-	//cout << "After pivot and update: " << endl << toString() << endl;
-}
-
-string SimplexWithBounds::toString() const
-{
-	string ret = "---\n" + m_state.toString();
-	for (auto const& [j, val]: m_assignments | ranges::views::enumerate)
-	{
-		ret += " " + m_state.variableNames[j];
-		if (m_basicVariables.count(j)) ret += " (b)";
-		ret += " = " + ::toString(val) + "\n";
-	}
-	ret += "------\n";
-	return ret;
-}
-
-pair<LPResult, variant<map<size_t, rational>, ReasonSet>> SolvingStateSimplifier::simplify()
-{
-	do
-	{
-		m_changed = false;
-		if (auto conflict = removeFixedVariables())
-			return {LPResult::Infeasible, move(*conflict)};
-		if (auto conflict = extractDirectConstraints())
-			return {LPResult::Infeasible, move(*conflict)};
-		// TODO we cannot do this anymore because it would
-		// mess up the variable numbering
-		// removeEmptyColumns();
-	}
-	while (m_changed);
-
-	return {LPResult::Unknown, move(m_fixedVariables)};
-}
-
-optional<ReasonSet> SolvingStateSimplifier::removeFixedVariables()
-{
-	for (auto const& [index, bounds]: m_state.bounds | ranges::views::enumerate)
-	{
-		if (!bounds.upper || (!bounds.lower && bounds.upper->numerator() > 0))
-			continue;
-		// Lower bound must be at least zero.
-		rational lower = max(rational{}, bounds.lower ? *bounds.lower : rational{});
-		rational upper = *bounds.upper;
-		if (upper < lower)
-			// Infeasible.
-			return bounds.lowerReasons + bounds.upperReasons;
-		if (upper != lower)
-			continue;
-		set<size_t> reasons = bounds.lowerReasons + bounds.upperReasons;
-		m_fixedVariables[index] = lower;
-		//cout << "Fixed " << m_state.variableNames.at(index) << " to " << ::toString(lower) << endl;
-		m_state.bounds[index] = {};
-		m_changed = true;
-
-		// substitute variable
-		// -> add the bounds to the literals for the conflict
-		// (maybe only if one of these constraints is used)
-		for (Constraint& constraint: m_state.constraints)
-			if (constraint.data[index])
-			{
-				constraint.data[0] -= constraint.data[index] * lower;
-				constraint.data[index] = 0;
-				constraint.reasons += reasons;
-			}
-	}
-
-	return nullopt;
-}
-
-optional<ReasonSet> SolvingStateSimplifier::extractDirectConstraints()
-{
-	vector<bool> constraintsToRemove(m_state.constraints.size(), false);
-	bool needsRemoval = false;
-	for (auto const& [index, constraint]: m_state.constraints | ranges::views::enumerate)
-	{
-		auto nonzeroCoefficients = constraint.data.enumerateTail() | ranges::views::filter(
-			[](std::pair<size_t, rational> const& _x) { return !!_x.second; }
-		);
-		// TODO we can exit early on in the loop above since we only care about zero, one or more than one nonzero entries.
-		// TODO could also use iterators and exit if we can advance it twice.
-		auto numNonzero = ranges::distance(nonzeroCoefficients);
-		if (numNonzero > 1)
-			continue;
-		constraintsToRemove[index] = true;
-		needsRemoval = true;
-		if (numNonzero == 0)
-		{
-			// 0 <= b or 0 = b
-			if (
-				constraint.data.front().numerator() < 0 ||
-				(constraint.equality && constraint.data.front())
-			)
-				return constraint.reasons;
-		}
-		else
-		{
-			auto&& [varIndex, factor] = nonzeroCoefficients.front();
-			// a * x <= b
-			rational bound = constraint.data[0] / factor;
-			if (
-				(factor >= 0 || constraint.equality) &&
-				(!m_state.bounds[varIndex].upper || bound < m_state.bounds[varIndex].upper)
-			)
-			{
-				m_state.bounds[varIndex].upper = bound;
-				m_state.bounds[varIndex].upperReasons = constraint.reasons;
-			}
-			if (
-				(factor <= 0 || constraint.equality) &&
-				bound >= 0 &&
-				(!m_state.bounds[varIndex].lower || bound > m_state.bounds[varIndex].lower)
-			)
-			{
-				m_state.bounds[varIndex].lower = bound;
-				m_state.bounds[varIndex].lowerReasons = constraint.reasons;
-			}
-		}
-	}
-	if (needsRemoval)
-	{
-		m_changed = true;
-		eraseIndices(m_state.constraints, constraintsToRemove);
-	}
-	return nullopt;
-}
-
-void SolvingStateSimplifier::removeEmptyColumns()
-{
-	vector<bool> variablesSeen(m_state.bounds.size(), false);
-	for (auto const& constraint: m_state.constraints)
-	{
-		for (auto&& [index, factor]: constraint.data.enumerateTail())
-			if (factor)
-				variablesSeen[index] = true;
-	}
-
-	vector<bool> variablesToRemove(variablesSeen.size(), false);
-	bool needsRemoval = false;
-	for (auto&& [i, seen]: variablesSeen | ranges::views::enumerate | ranges::views::tail)
-		if (!seen)
-		{
-			variablesToRemove[i] = true;
-			needsRemoval = true;
-			SolvingState::Bounds const& bounds = m_state.bounds.at(i);
-			// TODO actually it is unbounded if m_state.bounds.at(i).upper is nullopt.
-			if (bounds.lower || bounds.upper)
-			{
-				solAssert(!bounds.upper || bounds.upper >= 0);
-				if (bounds.lower && bounds.upper)
-					solAssert(*bounds.lower <= *bounds.upper);
-				m_fixedVariables[i] =
-					bounds.upper ?
-					*bounds.upper :
-					*bounds.lower;
-			}
-		}
-	if (needsRemoval)
-	{
-		m_changed = true;
-		removeColumns(m_state, variablesToRemove);
-	}
-}
-
-ProblemSplitter::ProblemSplitter(const SolvingState& _state):
-	m_state(_state),
-	m_column(1),
-	m_seenColumns(std::vector<bool>(m_state.variableNames.size(), false))
-{
-	while (m_column < m_state.variableNames.size() && nonZeroEntriesInColumn(_state, m_column).empty())
-		m_column++;
-}
-
-pair<vector<bool>, vector<bool>> ProblemSplitter::next()
-{
-	vector<bool> includedColumns;
-	vector<bool> includedRows;
-	tie(includedColumns, includedRows) = connectedComponent(m_state, m_column);
-
-	// Update state.
-	m_seenColumns |= includedColumns;
-	++m_column;
-	while (m_column < m_state.variableNames.size() && (
-		m_seenColumns[m_column] ||
-		nonZeroEntriesInColumn(m_state, m_column).empty()
-	))
-		++m_column;
-
-	return {move(includedColumns), move(includedRows)};
-/*
-	SolvingState splitOff;
-
-	splitOff.variableNames.emplace_back();
-	splitOff.bounds.emplace_back();
-
-	for (auto&& [i, included]: includedColumns | ranges::views::enumerate | ranges::views::tail)
-		if (included)
-		{
-			splitOff.variableNames.emplace_back(move(m_state.variableNames[i]));
-			splitOff.bounds.emplace_back(move(m_state.bounds[i]));
-		}
-
-	for (auto&& [i, included]: includedRows | ranges::views::enumerate)
-		if (included)
-		{
-			// Use const& on purpose because moving from the state can lead
-			// to undefined behaviour for connectedComponent
-			Constraint const& constraint = m_state.constraints[i];
-			Constraint splitRow{{}, constraint.equality, constraint.reasons};
-			for (size_t j = 0; j < constraint.data.size(); j++)
-				if (j == 0 || includedColumns[j])
-					splitRow.data.push_back(constraint.data[j]);
-			splitOff.constraints.push_back(move(splitRow));
-		}
-
-	return {includedColumns, splitOff};
-	*/
-}
-
-LPSolver::LPSolver(bool)
-{
-}
-
-LPResult LPSolver::setState(SolvingState _state)
-{
-	//cout << "Set state:\n" << _state.toString() << endl;
-	m_state = make_shared<SolvingState>(move(_state));
-	m_subProblems.clear();
-	m_subProblemsPerVariable.resize(m_state->variableNames.size(), static_cast<size_t>(-1));
-	m_subProblemsPerConstraint.resize(m_state->constraints.size(), static_cast<size_t>(-1));
-
-	normalizeRowLengths(*m_state);
-	auto&& [result, modelOrReasonSet] = SolvingStateSimplifier(*m_state).simplify();
-	if (result == LPResult::Infeasible)
-		return result;
-
-	// We do not need to store reasons because at this point, we do not have any reasons yet.
-	// We can add this and just need to store the reasons together with the variables.
-	m_fixedVariables = make_shared<std::map<size_t, rational>>(std::get<std::map<size_t, rational>>(modelOrReasonSet));
-
-	//cout << "Splitting..." << endl;
-	ProblemSplitter splitter(*m_state);
-	while (splitter)
-	{
-		auto&& [variables, constraints] = splitter.next();
-		SubProblem& problem = *m_subProblems.emplace_back(make_shared<SubProblem>());
-		solAssert(problem.dirty);
-		for (auto&& [i, included]: variables | ranges::views::enumerate)
-			if (included)
-			{
-				m_subProblemsPerVariable[i] = m_subProblems.size() - 1;
-				problem.variables.emplace(i);
-			}
-		for (auto&& [i, included]: constraints | ranges::views::enumerate)
-			if (included)
-				m_subProblemsPerConstraint[i] = m_subProblems.size() - 1;
-		//cout << "Adding new sub problem with " << m_subProblems.back()->variables.size() << " vars and " << m_subProblems.back()->constraints.size() << " constraints\n" << endl;
-		//cout << "-------------------- Split out" << endl;
-		//cout << m_subProblems.back()->state.toString() << endl;
-	}
-	//cout << "Done splitting." << endl;
-	return LPResult::Unknown;
-}
-
-void LPSolver::addConstraint(Constraint _constraint)
-{
-	//cout << "Adding constraint " << endl;
-	set<size_t> touchedProblems;
-	for (auto const& [index, entry]: _constraint.data.enumerateTail())
-		if (entry)
-		{
-			if (m_fixedVariables->count(index))
-			{
-				// This can directly lead to a conflict. We will check it later during the
-				// simplify run on the split problems.
-				_constraint.data[0] -= _constraint.data[index] * m_fixedVariables->at(index);
-				_constraint.data[index] = {};
-			}
-			else if (m_subProblemsPerVariable[index] != static_cast<size_t>(-1))
-				touchedProblems.emplace(m_subProblemsPerVariable[index]);
-		}
-	if (touchedProblems.empty())
-	{
-		//cout << "Creating new sub problem." << endl;
-		// TODO we could find an empty spot for the pointer.
-		m_subProblems.emplace_back(make_shared<SubProblem>());
-		solAssert(m_subProblems.back()->dirty);
-		touchedProblems.emplace(m_subProblems.size() - 1);
-	}
-	for (size_t problemToErase: touchedProblems | ranges::views::tail | ranges::views::reverse)
-		combineSubProblems(*touchedProblems.begin(), problemToErase);
-	addConstraintToSubProblem(*touchedProblems.begin(), _constraint);
-//	cout << "subproblems: " << (
-//		m_subProblems | ranges::views::transform([&](auto&& p) { return !!p; })).size() <<
-//		" (total: " << m_subProblems.size() << endl;
-	//cout << "done" << endl;
-}
-
-pair<LPResult, variant<Model, ReasonSet>> LPSolver::check()
-{
-	//cout << "Checking" << endl;
-	for (auto&& [index, problem]: m_subProblems | ranges::views::enumerate)
-	{
-		if (!problem)
-			continue;
-		if (!problem->dirty)
-		{
-			solAssert(problem->result != LPResult::Infeasible);
-			continue;
-		}
-
-		//cout << "Updating sub problem" << endl;
-
-		if (!problem->simplex)
-		{
-			auto&& [state, varMapping] = stateFromSubProblem(index);
-			problem->varMapping = move(varMapping);
-
-			// TODO if we do not use simplify, but just iteratively add constraints and bounsd
-			// to the simplex, we can detect "direct constraints" while adding (and compressing).
-			// and we might also be able to remove fixed variables.
-
-			// The simplify run is important because it detects conflicts
-			// due to fixed variables.
-			auto&& [result, modelOrReasonSet] = SolvingStateSimplifier{state}.simplify();
-			if (result == LPResult::Infeasible)
-			{
-				problem->result = LPResult::Infeasible;
-				problem->model = {};
-				problem->dirty = false;
-				// TODO we could use the improved reason set above.
-				return {LPResult::Infeasible, reasonSetForSubProblem(*problem)};
-			}
-			for (auto&& [fixedVar, value]: get<map<size_t, rational>>(modelOrReasonSet))
-			{
-				state.bounds.at(fixedVar).lower = value;
-				state.bounds.at(fixedVar).upper = value;
-			}
-
-			normalizeRowLengths(state);
-			problem->simplex = SimplexWithBounds{state};
-		}
-
-		// TODO we can just set the fixed variables to fixed bounds.
-		// Then adding new constraints using those is still fine.
-
-		//cout << state.toString() << endl;
-
-		// This is the new algorithm that uses bounds directly.
-		// TODO This new algorithm also allows us to lift the restriction
-		// that variables need to be non-negative.
-				/*
-		bool useSimplexWithBounds = true;
-		if (useSimplexWithBounds)
-		{*/
-			optional<LPResult> result;
-/*			if (m_cache)
-			{
-				auto it = m_cache->find(state);
-				if (it != m_cache->end())
-				{
-					//cout << "Cache hit" << endl;
-					result = it->second;
-				}
-			}*/
-			if (!result)
-			{
-				*result = problem->simplex->check();
-				// TODO we should even keep the updated tableau in the subproblem
-/*				if (m_cache)
-				{
-					(*m_cache)[state] = *result;
-					//cout << "Cache size " << m_cache->size() << endl;
-				}*/
-			}
-			problem->dirty = false;
-			problem->result = *result;
-			if (problem->result == LPResult::Infeasible)
-				return {LPResult::Infeasible, reasonSetForSubProblem(*problem)};
-		/*}
-		else if (auto conflict = boundsToConstraints(state))
-		{
-			problem->result = LPResult::Infeasible;
-			problem->model = {};
-			problem->dirty = false;
-			return {LPResult::Infeasible, reasonSetForSubProblem(*problem)};
-		}
-		else if (state.constraints.empty())
-		{
-			problem->result = LPResult::Feasible;
-			problem->model = {};
-			problem->dirty = false;
-		}
-		else
-		{
-			optional<LPResult> result;
-			if (m_cache)
-			{
-				auto it = m_cache->find(state);
-				if (it != m_cache->end())
-				{
-					//cout << "Cache hit" << endl;
-					result = it->second;
-				}
-			}
-			if (!result)
-			{
-				LinearExpression objectives;
-				objectives.resize(1);
-				objectives.resize(state.variableNames.size(), rational(bigint(1)));
-				// TODO the model relies on the variable numbering.
-				result = LPResult::Unknown;
-				tie(*result, problem->model) = simplex(state.constraints, move(objectives));
-				if (m_cache)
-				{
-					(*m_cache)[state] = *result;
-					//cout << "Cache size " << m_cache->size() << endl;
-				}
-			}
-			problem->dirty = false;
-			problem->result = *result;
-			if (problem->result == LPResult::Infeasible)
-				return {LPResult::Infeasible, reasonSetForSubProblem(*problem)};
-		}*/
-	}
-
-	return {LPResult::Unknown, Model{}};
-}
-
-void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom)
-{
-	//cout << "Combining " << _combineInto << " <- " << _combineFrom << endl;
-	// TOOD creating a copy and setting dirty is on operation.
-	m_subProblems[_combineInto] = make_shared<SubProblem>(*m_subProblems[_combineInto]);
-	m_subProblems[_combineInto]->dirty = true;
-	m_subProblems[_combineInto]->simplex = nullopt;
-	m_subProblems[_combineInto]->varMapping = {};
-
-	for (size_t& item: m_subProblemsPerVariable)
-		if (item == _combineFrom)
-			item = _combineInto;
-	for (size_t& item: m_subProblemsPerConstraint)
-		if (item == _combineFrom)
-			item = _combineInto;
-	m_subProblems[_combineInto]->variables += move(m_subProblems[_combineFrom]->variables);
-
-	m_subProblems[_combineFrom].reset();
-}
-
-void LPSolver::addConstraintToSubProblem(size_t _subProblem, Constraint _constraint)
-{
-	// TODO This is called in a loop - maybe we shouldn't copy every time.
-	m_subProblems[_subProblem] = make_shared<SubProblem>(*m_subProblems[_subProblem]);
-	SubProblem& problem = *m_subProblems[_subProblem];
-	problem.dirty = true;
-	if (problem.simplex)
-	{
-		Constraint compressedConstraint(_constraint);
-		compressedConstraint.data.resize(1, _constraint.data[0]);
-		for (auto const& [index, entry]: _constraint.data.enumerateTail())
-			if (entry)
-			{
-				size_t cIndex = 0;
-				if (problem.varMapping.count(index))
-				{
-					cIndex = problem.varMapping.at(index);
-					solAssert(m_subProblemsPerVariable.at(index) == _subProblem);
-				}
-				else
-				{
-					cIndex = problem.simplex->addVariable(m_state->variableNames.at(index));
-					problem.varMapping[index] = cIndex;
-					problem.variables.emplace(index);
-					solAssert(m_subProblemsPerVariable.at(index) == static_cast<size_t>(-1));
-					m_subProblemsPerVariable[index] = _subProblem;
-				}
-				solAssert(cIndex >= compressedConstraint.data.size());
-				compressedConstraint.data.resize(cIndex + 1);
-				compressedConstraint.data[cIndex] = entry;
-			}
-		problem.simplex->addConstraint(move(compressedConstraint));
-	}
-	else
-	{
-		problem.simplex = nullopt;
-		problem.varMapping = {};
-
-		for (auto const& [index, entry]: _constraint.data.enumerateTail())
-			if (entry)
-			{
-				solAssert(m_subProblemsPerVariable[index] == static_cast<size_t>(-1) || m_subProblemsPerVariable[index] == _subProblem);
-				m_subProblemsPerVariable[index] = _subProblem;
-				problem.variables.emplace(index);
-			}
-	}
-
-	problem.removableConstraints.emplace_back(move(_constraint));
-}
-
-pair<SolvingState, map<size_t, size_t>> LPSolver::stateFromSubProblem(size_t _index) const
-{
-	SolvingState split;
-	map<size_t, size_t> varMapping;
-
-	split.variableNames.emplace_back();
-	split.bounds.emplace_back();
-
-	SubProblem const& problem = *m_subProblems[_index];
-	for (size_t varIndex: problem.variables)
-	{
-		varMapping[varIndex] = split.variableNames.size();
-		split.variableNames.emplace_back(m_state->variableNames[varIndex]);
-		split.bounds.emplace_back(m_state->bounds[varIndex]);
-	}
-	for (auto&& item: m_subProblemsPerConstraint | ranges::views::enumerate)
-		if (item.second == _index)
-		{
-			Constraint const& constraint = m_state->constraints[item.first];
-			Constraint splitRow{{}, constraint.equality, constraint.reasons};
-			splitRow.data.push_back(constraint.data.get(0));
-			for (size_t varIndex: problem.variables)
-				splitRow.data.push_back(constraint.data.get(varIndex));
-			split.constraints.push_back(move(splitRow));
-		}
-
-	for (Constraint const& constraint: m_subProblems[_index]->removableConstraints)
-	{
-		Constraint splitRow{{}, constraint.equality, constraint.reasons};
-		splitRow.data.push_back(constraint.data.get(0));
-		for (size_t varIndex: problem.variables)
-			splitRow.data.push_back(constraint.data.get(varIndex));
-		split.constraints.push_back(move(splitRow));
-	}
-
-	normalizeRowLengths(split);
-
-	return {split, varMapping};
-}
-
-ReasonSet LPSolver::reasonSetForSubProblem(LPSolver::SubProblem const& _subProblem)
-{
-	ReasonSet reasons;
-	for (Constraint const& c: _subProblem.removableConstraints)
-		reasons += c.reasons;
-	return reasons;
 }
diff --git a/libsolutil/LP.h b/libsolutil/LP.h
index 56daeb2fe..c1ea8734d 100644
--- a/libsolutil/LP.h
+++ b/libsolutil/LP.h
@@ -42,8 +42,6 @@ struct Constraint
 {
 	LinearExpression data;
 	bool equality = false;
-	/// Set of literals the conjunction of which implies this constraint.
-	std::set<size_t> reasons = {};
 
 	bool operator<(Constraint const& _other) const;
 	bool operator==(Constraint const& _other) const;
@@ -171,156 +169,92 @@ enum class LPResult
 	Infeasible ///< System does not have any solution.
 };
 
-
-class SimplexWithBounds
-{
-public:
-	explicit SimplexWithBounds(SolvingState _state);
-	LPResult check();
-
-	size_t addVariable(std::string _name);
-	void addConstraint(Constraint _constraint);
-
-	std::string toString() const;
-private:
-	/// Set value of non-basic variable.
-	void update(size_t _var, rational const& _value);
-	/// @returns the index of the first basic variable violating its bounds.
-	std::optional<size_t> firstConflictingBasicVariable() const;
-	std::optional<size_t> firstReplacementVar(size_t _basicVarToReplace, bool _increasing) const;
-
-	void pivot(size_t _old, size_t _new);
-	void pivotAndUpdate(size_t _oldBasicVar, rational const& _newValue, size_t _newBasicVar);
-
-	SolvingState m_state;
-	std::vector<rational> m_assignments;
-	/// Variable index to row it controls.
-	std::map<size_t, size_t> m_basicVariables;
-};
-
-
 /**
- * Applies several strategies to simplify a given solving state.
- * During these simplifications, it can sometimes already be determined if the
- * state is feasible or not.
- * Since some variables can be fixed to specific values, it returns a
- * (partial) model.
- *
- * - Constraints with exactly one nonzero coefficient represent "a x <= b"
- *   and thus are turned into bounds.
- * - Constraints with zero nonzero coefficients are constant relations.
- *   If such a relation is false, answer "infeasible", otherwise remove the constraint.
- * - Empty columns can be removed.
- * - Variables with matching bounds can be removed from the problem by substitution.
- *
- * Holds a reference to the solving state that is modified during operation.
- */
-class SolvingStateSimplifier
-{
-public:
-	SolvingStateSimplifier(SolvingState& _state):
-		m_state(_state) {}
-
-	std::pair<LPResult, std::variant<std::map<size_t, rational>, ReasonSet>> simplify();
-
-private:
-	/// Remove variables that have equal lower and upper bound.
-	/// @returns reason / set of conflicting clauses if infeasible.
-	std::optional<ReasonSet> removeFixedVariables();
-
-	/// Removes constraints of the form 0 <= b or 0 == b (no variables) and
-	/// turns constraints of the form a * x <= b (one variable) into bounds.
-	/// @returns reason / set of conflicting clauses if infeasible.
-	std::optional<ReasonSet> extractDirectConstraints();
-
-	/// Removes all-zeros columns.
-	void removeEmptyColumns();
-
-	/// Set to true by the strategies if they performed some changes.
-	bool m_changed = false;
-
-	SolvingState& m_state;
-	std::map<size_t, rational> m_fixedVariables;
-};
-
-/**
- * Splits a given linear program into multiple linear programs with disjoint sets of variables.
- * The initial program is feasible if and only if all sub-programs are feasible.
- */
-class ProblemSplitter
-{
-public:
-	explicit ProblemSplitter(SolvingState const& _state);
-
-	/// @returns true if there are still sub-problems to split out.
-	operator bool() const { return m_column < m_state.variableNames.size(); }
-
-	/// @returns the next sub-problem.
-	std::pair<std::vector<bool>, std::vector<bool>> next();
-
-private:
-	SolvingState const& m_state;
-	/// Next column to start the search for a connected component.
-	size_t m_column = 1;
-	/// The columns we have already split out.
-	std::vector<bool> m_seenColumns;
-};
-
-/**
- * LP solver for rational problems.
+ * LP solver for rational problems, based on "A Fast Linear-Arithmetic Solver for DPLL(T)*"
+ * by Dutertre and Moura.
  *
  * Does not solve integer problems!
  *
- * Tries to split a given problem into sub-problems and utilizes a cache to quickly solve
- * similar problems.
+ * Tries to split incoming bounds and constraints into unrelated sub-problems.
+ * Maintains lower/upper bounds for all variables.
+ * Adds one slack variable per constraint and stores all constraints as "= 0" equations.
+ * Splits variables into basic and non-basic. For each row there is exactly one
+ * basic variable that has a factor of -1.
+ * The equations are satisfied at all times and non-basic variables are always within their bounds.
+ * Non-basic variables might violate their bounds.
+ * It attempts to resolve these violations in turn, swapping a basic variables with a non-basic
+ * variables that can still move in the required direction.
  *
- * Can be used in a mode where it does not support returning models. In that case, the
- * cache is more efficient.
+ * It is perfectly fine to add new bounds, variable or constraints after a call to "check".
+ * The solver can be copied at low cost and it uses a "copy on write" mechanism for the sub-problems.
  */
 class LPSolver
 {
 public:
-	explicit LPSolver(bool _supportModels = true);
-	explicit LPSolver(std::unordered_map<SolvingState, LPResult>* _cache):
-		m_cache(_cache) {}
+	void addConstraint(Constraint const& _constraint, std::optional<size_t> _reason = std::nullopt);
+	void setVariableName(size_t _variable, std::string _name);
+	void addLowerBound(size_t _variable, rational _bound);
+	void addUpperBound(size_t _variable, rational _bound);
 
-
-	LPResult setState(SolvingState _state);
-	void addConstraint(Constraint _constraint);
 	std::pair<LPResult, std::variant<Model, ReasonSet>> check();
 
-private:
-	void combineSubProblems(size_t _combineInto, size_t _combineFrom);
-	void addConstraintToSubProblem(size_t _subProblem, Constraint _constraint);
-	void updateSubProblems();
+	std::string toString() const;
 
-	/// Ground state for CDCL. This is shared by copies of the solver.
-	/// Only ``setState`` changes the state. Copies will only use
-	/// ``addConstraint`` which does not change m_state.
-	std::shared_ptr<SolvingState> m_state;
+private:
+	struct Bounds
+	{
+		std::optional<rational> lower;
+		std::optional<rational> upper;
+	};
+	struct Variable
+	{
+		std::string name = {};
+		rational value = 0;
+		Bounds bounds = {};
+	};
 	struct SubProblem
 	{
-		// TODO now we could actually put the constraints here again.
-		std::vector<Constraint> removableConstraints;
-		bool dirty = true;
-		LPResult result = LPResult::Unknown;
-		std::vector<boost::rational<bigint>> model = {};
-		std::set<size_t> variables = {};
-		std::optional<SimplexWithBounds> simplex = std::nullopt;
+		/// Set to true on "check". Needs a copy for adding a constraint or bound if set to true.
+		bool sealed = false;
+		std::optional<LPResult> result = std::nullopt;
+		std::vector<LinearExpression> factors;
+		std::vector<Variable> variables;
+		/// Variable index to constraint it controls.
+		std::map<size_t, size_t> basicVariables;
+		/// Maps outer indices to inner indices.
 		std::map<size_t, size_t> varMapping = {};
+		std::set<size_t> reasons;
+
+		LPResult check();
+		std::string toString() const;
+	private:
+		/// Set value of non-basic variable.
+		void update(size_t _varIndex, rational const& _value);
+		/// @returns the index of the first basic variable violating its bounds.
+		std::optional<size_t> firstConflictingBasicVariable() const;
+		std::optional<size_t> firstReplacementVar(size_t _basicVarToReplace, bool _increasing) const;
+
+		void pivot(size_t _old, size_t _new);
+		void pivotAndUpdate(size_t _oldBasicVar, rational const& _newValue, size_t _newBasicVar);
 	};
 
-	std::pair<SolvingState, std::map<size_t, size_t>> stateFromSubProblem(size_t _index) const;
-	ReasonSet reasonSetForSubProblem(SubProblem const& _subProblem);
 
-	std::shared_ptr<std::map<size_t, rational>> m_fixedVariables;
+	SubProblem& unseal(size_t _problem);
+	/// Unseals the problem for the given variable or creates a new one.
+	SubProblem& unsealForVariable(size_t _outerIndex);
+	void combineSubProblems(size_t _combineInto, size_t _combineFrom);
+	void addConstraintToSubProblem(size_t _subProblem, Constraint const& _constraint, std::optional<size_t> _reason);
+	void addOuterVariableToSubProblem(size_t _subProblem, size_t _outerIndex);
+	size_t addNewVariableToSubProblem(size_t _subProblem);
+
+	std::map<std::string, rational> model() const;
+
 	/// These use "copy on write".
 	std::vector<std::shared_ptr<SubProblem>> m_subProblems;
-	std::vector<size_t> m_subProblemsPerVariable;
-	std::vector<size_t> m_subProblemsPerConstraint;
-	/// TODO also store the first infeasible subproblem?
-	/// TODO still retain the cache?
-	std::unordered_map<SolvingState, LPResult>* m_cache = nullptr;
+	/// Maps outer indices to sub problems.
+	std::map<size_t, size_t> m_subProblemsPerVariable;
+	/// Counter to enable unique names for the slack variables.
+	size_t m_slackVariableCounter = 0;
 
 };
 
diff --git a/test/libsolutil/LP.cpp b/test/libsolutil/LP.cpp
index e9e05de1a..331190479 100644
--- a/test/libsolutil/LP.cpp
+++ b/test/libsolutil/LP.cpp
@@ -38,7 +38,6 @@ class LPTestFramework
 public:
 	LPTestFramework()
 	{
-		m_solvingState.variableNames.emplace_back("");
 	}
 
 	LinearExpression constant(rational _value)
@@ -52,11 +51,11 @@ public:
 	}
 
 	/// Adds the constraint "_lhs <= _rhs".
-	void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs, set<size_t> _reason = {})
+	void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs, optional<size_t> _reason = {})
 	{
 		_lhs -= _rhs;
 		_lhs[0] = -_lhs[0];
-		m_solvingState.constraints.push_back({move(_lhs), false, move(_reason)});
+		m_solver.addConstraint({move(_lhs), false}, move(_reason));
 	}
 
 	void addLEConstraint(LinearExpression _lhs, rational _rhs)
@@ -65,47 +64,43 @@ public:
 	}
 
 	/// Adds the constraint "_lhs = _rhs".
-	void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs, set<size_t> _reason = {})
+	void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs, optional<size_t> _reason = {})
 	{
 		_lhs -= _rhs;
 		_lhs[0] = -_lhs[0];
-		m_solvingState.constraints.push_back({move(_lhs), true, move(_reason)});
+		m_solver.addConstraint({move(_lhs), true}, move(_reason));
 	}
 
-	void addLowerBound(string _variable, rational _value, set<size_t> _reason = {})
+	void addLowerBound(string _variable, rational _value)
 	{
-		size_t index = variableIndex(_variable);
-		if (index >= m_solvingState.bounds.size())
-			m_solvingState.bounds.resize(index + 1);
-		m_solvingState.bounds.at(index).lower = _value;
-		m_solvingState.bounds.at(index).lowerReasons = move(_reason);
+		m_solver.addLowerBound(variableIndex(_variable), _value);
 	}
 
-	void addUpperBound(string _variable, rational _value, set<size_t> _reason = {})
+	void addUpperBound(string _variable, rational _value)
 	{
-		size_t index = variableIndex(_variable);
-		if (index >= m_solvingState.bounds.size())
-			m_solvingState.bounds.resize(index + 1);
-		m_solvingState.bounds.at(index).upper = _value;
-		m_solvingState.bounds.at(index).upperReasons = move(_reason);
+		m_solver.addUpperBound(variableIndex(_variable), _value);
 	}
 
 	void feasible(vector<pair<string, rational>> const& _solution)
 	{
-		m_solver.setState(m_solvingState);
 		auto [result, modelOrReasonSet] = m_solver.check();
 		BOOST_REQUIRE(result == LPResult::Feasible);
 		Model model = get<Model>(modelOrReasonSet);
 		for (auto const& [var, value]: _solution)
+		{
+			BOOST_CHECK_MESSAGE(
+				model.count(var),
+				"Variable " + var + " not found in model."
+			);
 			BOOST_CHECK_MESSAGE(
 				value == model.at(var),
 				var + " = "s + ::toString(model.at(var)) + " (expected " + ::toString(value) + ")"
 			);
+		}
 	}
 
 	void infeasible(set<size_t> _reason = {})
 	{
-		m_solver.setState(m_solvingState);
 		auto [result, modelOrReason] = m_solver.check();
 		BOOST_REQUIRE(result == LPResult::Infeasible);
 		ReasonSet suppliedReason = get<ReasonSet>(modelOrReason);
@@ -115,44 +110,54 @@ public:
 protected:
 	size_t variableIndex(string const& _name)
 	{
-		if (m_solvingState.variableNames.empty())
-			m_solvingState.variableNames.emplace_back("");
-		auto index = findOffset(m_solvingState.variableNames, _name);
-		if (!index)
+		if (!m_variableIndices.count(_name))
 		{
-			index = m_solvingState.variableNames.size();
-			m_solvingState.variableNames.emplace_back(_name);
+			size_t index = 1 + m_variableIndices.size();
+			m_solver.setVariableName(index, _name);
+			m_variableIndices[_name] = index;
 		}
-		return *index;
+		return m_variableIndices.at(_name);
 	}
 
 	LPSolver m_solver;
-	SolvingState m_solvingState;
+	map<string, size_t> m_variableIndices;
 };
 
 
 BOOST_FIXTURE_TEST_SUITE(LP, LPTestFramework, *boost::unit_test::label("nooptions"))
 
+BOOST_AUTO_TEST_CASE(empty)
+{
+	feasible({});
+}
+
 BOOST_AUTO_TEST_CASE(basic)
 {
 	auto x = variable("x");
 	addLEConstraint(2 * x, 10);
-	feasible({{"x", 5}});
+	feasible({{"x", 0}});
 }
 
 BOOST_AUTO_TEST_CASE(not_linear_independent)
 {
 	addLEConstraint(2 * variable("x"), 10);
 	addLEConstraint(4 * variable("x"), 20);
+	feasible({{"x", 0}});
+}
+
+BOOST_AUTO_TEST_CASE(not_linear_independent_eq)
+{
+	addLEConstraint(2 * variable("x"), 10);
+	addEQConstraint(4 * variable("x"), constant(20));
 	feasible({{"x", 5}});
 }
 
 BOOST_AUTO_TEST_CASE(two_vars)
 {
 	addLEConstraint(variable("y"), 3);
-	addLEConstraint(variable("x"), 10);
+	addLowerBound("x", 10);
 	addLEConstraint(variable("x") + variable("y"), 4);
-	feasible({{"x", 1}, {"y", 3}});
+	feasible({{"x", 10}, {"y", -6}});
 }
 
 BOOST_AUTO_TEST_CASE(one_le_the_other)
@@ -167,8 +172,8 @@ BOOST_AUTO_TEST_CASE(factors)
 	auto y = variable("y");
 	addLEConstraint(2 * y, 3);
 	addLEConstraint(16 * x, 10);
-	addLEConstraint(x + y, 4);
-	feasible({{"x", rational(5) / 8}, {"y", rational(3) / 2}});
+	addLEConstraint(constant(2), x + y);
+	feasible({{"x", rational(5) / 8}, {"y", rational(11) / 8}});
 }
 
 BOOST_AUTO_TEST_CASE(cache)
@@ -178,8 +183,8 @@ BOOST_AUTO_TEST_CASE(cache)
 	// that it results in the same value.
 	auto x = variable("x");
 	auto y = variable("y");
-	addLEConstraint(2 * y, 3);
-	addLEConstraint(2 * x, 3);
+	addLEConstraint(constant(3), 2 * y);
+	addLEConstraint(constant(3), 2 * x);
 	feasible({{"x", rational(3) / 2}, {"y", rational(3) / 2}});
 	feasible({{"x", rational(3) / 2}, {"y", rational(3) / 2}});
 }
@@ -187,19 +192,25 @@ BOOST_AUTO_TEST_CASE(cache)
 BOOST_AUTO_TEST_CASE(bounds)
 {
 	addUpperBound("x", 200);
-	feasible({{"x", 200}});
+	feasible({{"x", 0}});
 
 	addLEConstraint(variable("x"), 100);
-	feasible({{"x", 100}});
+	feasible({{"x", 0}});
 
 	addLEConstraint(constant(5), variable("x"));
-	feasible({{"x", 100}});
+	feasible({{"x", 5}});
 
 	addLowerBound("x", 20);
-	feasible({{"x", 100}});
+	feasible({{"x", 20}});
 	addLowerBound("x", 25);
-	feasible({{"x", 100}});
+	feasible({{"x", 25}});
+	addLowerBound("x", 21);
+	feasible({{"x", 25}});
 
+	addUpperBound("x", 26);
+	feasible({{"x", 25}});
+	addUpperBound("x", 28);
+	feasible({{"x", 25}});
 	addUpperBound("x", 20);
 	infeasible();
 }
@@ -210,19 +221,20 @@ BOOST_AUTO_TEST_CASE(bounds2)
 	addUpperBound("x", 250);
 	addLowerBound("y", 2);
 	addUpperBound("y", 3);
-	feasible({{"x", 250}, {"y", 3}});
+	feasible({{"x", 200}, {"y", 2}});
 
 	addLEConstraint(variable("y"), variable("x"));
-	feasible({{"x", 250}, {"y", 3}});
+	feasible({{"x", 200}, {"y", 2}});
 
 	addEQConstraint(variable("y") + constant(231), variable("x"));
-	feasible({{"x", 234}, {"y", 3}});
-
+	cout << m_solver.toString() << endl;
+	feasible({{"x", 233}, {"y", 2}});
+/*
 	addEQConstraint(variable("y") + constant(10), variable("x") - variable("z"));
-	feasible({{"x", 234}, {"y", 3}});
+	feasible({{"x", 234}, {"y", 3}, {"z", 0}});
 
 	addEQConstraint(variable("z") + variable("x"), constant(2));
-	infeasible();
+	infeasible();*/
 }
 
 BOOST_AUTO_TEST_CASE(lower_bound)
@@ -230,7 +242,7 @@ BOOST_AUTO_TEST_CASE(lower_bound)
 	addLEConstraint(constant(1), variable("y"));
 	addLEConstraint(variable("x"), constant(10));
 	addLEConstraint(2 * variable("x") + variable("y"), 2);
-	feasible({{"x", 0}, {"y", 2}});
+	feasible({{"x", 0}, {"y", 1}});
 }
 
 BOOST_AUTO_TEST_CASE(check_infeasible)
@@ -252,11 +264,13 @@ BOOST_AUTO_TEST_CASE(unbounded2)
 	auto y = variable("y");
 	addLEConstraint(constant(2), x + y);
 	addLEConstraint(x, 10);
-	feasible({{"x", 10}, {"y", 0}});
+	feasible({{"x", 2}, {"y", 0}});
 }
 
 BOOST_AUTO_TEST_CASE(unbounded3)
 {
+	addLowerBound("y", 0);
+	addLowerBound("x", 0);
 	addLEConstraint(constant(0) - variable("x") - variable("y"), constant(10));
 	feasible({{"x", 0}, {"y", 0}});
 
@@ -277,7 +291,7 @@ BOOST_AUTO_TEST_CASE(equal)
 	auto y = variable("y");
 	addEQConstraint(x, y + constant(10));
 	addLEConstraint(x, 20);
-	feasible({{"x", 20}, {"y", 10}});
+	feasible({{"x", 10}, {"y", 0}});
 }
 
 
@@ -285,6 +299,7 @@ BOOST_AUTO_TEST_CASE(equal_constant)
 {
 	auto x = variable("x");
 	auto y = variable("y");
+	addLowerBound("x", 5);
 	addLEConstraint(x, y);
 	addEQConstraint(y, constant(5));
 	feasible({{"x", 5}, {"y", 5}});
@@ -295,7 +310,26 @@ BOOST_AUTO_TEST_CASE(all_equality)
 	auto x = variable("x");
 	auto y = variable("y");
 	addEQConstraint(-6 * x - 6 * y, constant(8));
-	infeasible();
+	feasible({{"x", -rational(4) / 3}, {"y", 0}});
+}
+
+BOOST_AUTO_TEST_CASE(negative_values)
+{
+	auto x = variable("x");
+	auto y = variable("y");
+	addEQConstraint(-2 * x, constant(8));
+	addLowerBound("y", -2);
+	addUpperBound("y", -1);
+	feasible({{"x", -4}, {"y", -1}});
+}
+
+BOOST_AUTO_TEST_CASE(basic_basic)
+{
+	auto x = variable("x");
+	addLEConstraint(x, constant(5));
+	feasible({{"x", 0}});
+	addEQConstraint(x, constant(5));
+	feasible({{"x", 5}});
 }
 
 BOOST_AUTO_TEST_CASE(linear_dependent)
@@ -303,40 +337,40 @@ BOOST_AUTO_TEST_CASE(linear_dependent)
 	auto x = variable("x");
 	auto y = variable("y");
 	auto z = variable("z");
-	addLEConstraint(x, 5);
-	addLEConstraint(2 * y, 10);
-	addLEConstraint(3 * z, 15);
+	addLEConstraint(x, -5);
+	addLEConstraint(2 * y, -10);
+	addLEConstraint(3 * z, -15);
 	// Here, they should be split into three independent problems.
-	feasible({{"x", 5}, {"y", 5}, {"z", 5}});
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
 	addLEConstraint((x + y) + z, 100);
-	feasible({{"x", 5}, {"y", 5}, {"z", 5}});
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
-	addLEConstraint((x + y) + z, 2);
-	feasible({{"x", 2}, {"y", 0}, {"z", 0}});
+	addLEConstraint((x + y) + z, -1);
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
-	addLEConstraint(constant(2), (x + y) + z);
-	feasible({{"x", 2}, {"y", 0}, {"z", 0}});
+	addLEConstraint(constant(-20), (x + y) + z);
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
-	addEQConstraint(constant(2), (x + y) + z);
-	feasible({{"x", 2}, {"y", 0}, {"z", 0}});
+	addEQConstraint(constant(-18), (x + y) + z);
+	feasible({{"x", -8}, {"y", -5}, {"z", -5}});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_simple)
 {
 	auto x = variable("x");
 	addLEConstraint(2 * x, constant(20), {0});
-	addLowerBound("x", 12, {1});
-	infeasible({0, 1});
+	addLowerBound("x", 12);
+	infeasible({0});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_bounds)
 {
 	auto x = variable("x");
-	addLowerBound("x", 12, {0});
-	addUpperBound("x", 11, {1});
-	addLEConstraint(x, constant(200), {2}); // unrelated
-	infeasible({0, 1});
+	addLowerBound("x", 12);
+	addUpperBound("x", 11);
+	addLEConstraint(variable("y"), constant(200), {2}); // unrelated
+	infeasible({});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_forwarded)
@@ -346,9 +380,9 @@ BOOST_AUTO_TEST_CASE(reasons_forwarded)
 	addEQConstraint(x, constant(2), {0});
 	addEQConstraint(x, y, {1});
 	feasible({});
-	addLEConstraint(x, constant(200), {5}); // unrelated
+	addLEConstraint(x, constant(200), {5});
 	addLEConstraint(y, constant(1), {2});
-	infeasible({0, 1, 2});
+	infeasible({0, 1, 2, 5});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_forwarded2)
@@ -359,9 +393,9 @@ BOOST_AUTO_TEST_CASE(reasons_forwarded2)
 	addEQConstraint(x, y, {1});
 	feasible({});
 	addEQConstraint(y, constant(3), {2});
-	addLEConstraint(x, constant(200), {6}); // unrelated
-	addLEConstraint(y, constant(202), {6}); // unrelated
-	infeasible({0, 1, 2});
+	addLEConstraint(x, constant(200), {6});
+	addLEConstraint(y, constant(202), {7});
+	infeasible({0, 1, 2, 6, 7});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_split)
@@ -377,6 +411,22 @@ BOOST_AUTO_TEST_CASE(reasons_split)
 	infeasible({0, 2});
 }
 
+BOOST_AUTO_TEST_CASE(reasons_joined)
+{
+	auto x = variable("x");
+	auto y = variable("y");
+	auto z = variable("z");
+	addLEConstraint(x + y, constant(2), {0});
+	feasible({});
+	addLEConstraint(constant(20), x + y, {2});
+	infeasible({0, 2});
+	addLEConstraint(z, constant(200), {3});
+	infeasible({0, 2});
+	addLEConstraint(x, z);
+	infeasible({0, 2, 3});
+}
+
+
 BOOST_AUTO_TEST_CASE(fuzzer2)
 {
 	/*
@@ -399,6 +449,15 @@ BOOST_AUTO_TEST_CASE(fuzzer2)
 	auto x6 = variable("x6");
 	auto x7 = variable("x7");
 	auto x8 = variable("x8");
+	addLowerBound("x0", 0);
+	addLowerBound("x1", 0);
+	addLowerBound("x2", 0);
+	addLowerBound("x3", 0);
+	addLowerBound("x4", 0);
+	addLowerBound("x5", 0);
+	addLowerBound("x6", 0);
+	addLowerBound("x7", 0);
+	addLowerBound("x8", 0);
 	addEQConstraint(90 * x4 + 9 * x5, constant(-1));
 	addLEConstraint(-76 * x2 + 74 * x5, constant(0));
 	addLEConstraint(31 * x3 - 71 * x8, constant(0));