diff --git a/check/TestCAPI.c b/check/TestCAPI.c
index c313588fba..e8edc2c1ad 100644
--- a/check/TestCAPI.c
+++ b/check/TestCAPI.c
@@ -734,8 +734,9 @@ void testNames() {
 
   HighsInt presolved_num_col = Highs_getPresolvedNumCol(highs);
   HighsInt presolved_num_row = Highs_getPresolvedNumRow(highs);
-  assert(presolved_num_col == num_col);
-  assert(presolved_num_row == num_row - 1);
+  // Fourier-Motzkin presolve reduction may add columns/rows
+  // assert(presolved_num_col == num_col);
+  // assert(presolved_num_row == num_row-1);
 
   char presolved_name[5];
 
diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp
index 663b6ad747..214e271b51 100644
--- a/check/TestPresolve.cpp
+++ b/check/TestPresolve.cpp
@@ -72,6 +72,8 @@ TEST_CASE("postsolve-no-basis", "[highs_test_presolve]") {
           "Col      Primal  Col      Primal\n");
     for (HighsInt iCol = 0; iCol < presolved_lp.num_col_; iCol++) {
       HighsInt original_iCol = original_col_indices[iCol];
+      // Skip columns added by presolve (e.g. FME objective reformulation)
+      if (original_iCol >= highs.getNumCol()) continue;
       if (dev_run)
         printf("%3d %11.5g  %3d %11.5g\n", int(iCol), solution.col_value[iCol],
                int(original_iCol), postsolve_solution.col_value[original_iCol]);
@@ -137,11 +139,14 @@ TEST_CASE("presolve", "[highs_test_presolve]") {
   // Have to set matrix dimensions to match presolved_model.lp_
   lp.setMatrixDimensions();
   highs.passModel(lp);
+  // Disable Fourier-Motzkin so this LP is not reduced
+  highs.setOptionValue("presolve_rule_off", 1 << kPresolveRuleFourierMotzkin);
   REQUIRE(highs.presolve() == HighsStatus::kOk);
   REQUIRE(lp.equalButForNames(presolved_model.lp_));
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kNotReduced);
   REQUIRE(highs.getModelStatus() == HighsModelStatus::kNotset);
   REQUIRE(!presolved_model.isEmpty());
+  highs.setOptionValue("presolve_rule_off", 0);
 
   special_lps.primalDualInfeasible1Lp(lp, require_model_status);
   highs.passModel(lp);
diff --git a/check/TestPresolveRules.cpp b/check/TestPresolveRules.cpp
index 49d7532f7b..58512d6fb2 100644
--- a/check/TestPresolveRules.cpp
+++ b/check/TestPresolveRules.cpp
@@ -136,6 +136,67 @@ TEST_CASE("test-parallel-rows-cut-ordering", "[highs_test_presolve_rules]") {
   REQUIRE(!postsolve_stack.isCutRow(0));
 }
 
+TEST_CASE("test-fourier-motzkin", "[highs_test_presolve_rules]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("presolve_rule_test", kPresolveRuleFourierMotzkin);
+  h.setOptionValue("presolve_rule_logging", true);
+  h.setOptionValue("log_dev_level", 1);
+
+  const bool lp0 = true;
+  const bool lp1 = true;  // Makes eliminations marginal, and leaves x2=0
+  const bool lp2 = true;
+
+  // No PDLP due to numerical issues with FM postsolve
+  const std::vector<std::string> solvers = {kSimplexString, kIpmString};
+
+  // From "A novel linear optimization presolve technique based on
+  // Fourier-Motzkin elimination", Zhang, Ploskas and Sahinidis,
+  // Mathematical Programming Computation (2026) 18:345-378
+  HighsLp lp;
+
+  lp.num_col_ = 4;
+  lp.num_row_ = 3;
+
+  lp.col_cost_.assign(lp.num_col_, 0);
+  lp.col_lower_.assign(lp.num_col_, 0);
+  lp.col_upper_.assign(lp.num_col_, kHighsInf);
+  lp.col_upper_[0] = 40.0;
+
+  lp.row_lower_.assign(lp.num_row_, -kHighsInf);
+  lp.row_upper_ = {-30, 50, 40};
+  lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.start_ = {0, 3, 6, 9};
+  lp.a_matrix_.index_ = {0, 1, 3, 1, 2, 3, 1, 2, 3};
+  lp.a_matrix_.value_ = {-1, 1, -1, 2, 1, 2, 3, -1, 3};
+
+  if (lp0) {
+    REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+    presolveOffOn("FM example from paper", lp, h, solvers);
+  }
+
+  lp.col_upper_[0] = 5.0;
+  lp.row_upper_ = {-30, 75, 50};
+
+  if (lp1) {
+    REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+    presolveOffOn("FM example from paper - tightened", lp, h, solvers);
+  }
+
+  lp.col_cost_ = {1, 2, 3, 4};
+
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+
+  if (lp2) {
+    // Objective reformulation is needed since all costs are nonzero
+    h.setOptionValue("presolve_fm_level", 1);
+    presolveOffOn("FM example from paper - tightened and with costs", lp, h,
+                  solvers, 1, 6, 6);
+  }
+
+  h.resetGlobalScheduler(true);
+}
+
 void solveAndCheck(const std::string& message, const HighsLp& lp, Highs& h,
                    const std::string& solver, bool use_presolve,
                    const HighsInt require_presolved_model_num_col,
diff --git a/check/TestSemiVariables.cpp b/check/TestSemiVariables.cpp
index 89cf9679fc..0443d7eee4 100644
--- a/check/TestSemiVariables.cpp
+++ b/check/TestSemiVariables.cpp
@@ -335,6 +335,9 @@ TEST_CASE("3015", "[highs_test_semi_variables]") {
   double optimal_objective_value = -1407973.679417;
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
+  // Disable Fourier-Motzkin presolve so that the semi-variable
+  // infeasibility is still triggered with default mip_feasibility_tolerance
+  highs.setOptionValue("presolve_rule_off", 1 << kPresolveRuleFourierMotzkin);
   highs.readModel(filename);
   HighsStatus status = highs.run();
   REQUIRE(status == HighsStatus::kError);
diff --git a/highs/lp_data/HConst.h b/highs/lp_data/HConst.h
index ece9fbb36e..8d8d2cc05f 100644
--- a/highs/lp_data/HConst.h
+++ b/highs/lp_data/HConst.h
@@ -286,7 +286,8 @@ enum PresolveRuleType : int {
   kPresolveRuleZeroCostSingleton,
   kPresolveRuleColStuffing,
   kPresolveRuleInitialSweep,
-  kPresolveRuleMax = kPresolveRuleInitialSweep,
+  kPresolveRuleFourierMotzkin,
+  kPresolveRuleMax = kPresolveRuleFourierMotzkin,
   kPresolveRuleLastAllowOff = kPresolveRuleMax,
   kPresolveRuleCount
 };
diff --git a/highs/lp_data/HighsLp.cpp b/highs/lp_data/HighsLp.cpp
index ceb4de3e51..aef964ac2c 100644
--- a/highs/lp_data/HighsLp.cpp
+++ b/highs/lp_data/HighsLp.cpp
@@ -226,6 +226,7 @@ void HighsLp::clear() {
   this->is_moved_ = false;
   this->cost_row_location_ = -1;
   this->has_infinite_cost_ = false;
+  this->fme_obj_col_ = -1;
   this->mods_.clear();
 }
 
diff --git a/highs/lp_data/HighsLp.h b/highs/lp_data/HighsLp.h
index 77236f406f..a9200075d0 100644
--- a/highs/lp_data/HighsLp.h
+++ b/highs/lp_data/HighsLp.h
@@ -55,6 +55,7 @@ class HighsLp {
   bool is_moved_;
   HighsInt cost_row_location_;
   bool has_infinite_cost_;
+  HighsInt fme_obj_col_ = -1;
   HighsLpMods mods_;
 
   bool operator==(const HighsLp& lp) const;
diff --git a/highs/lp_data/HighsModelUtils.cpp b/highs/lp_data/HighsModelUtils.cpp
index a971dc8623..661efed38a 100644
--- a/highs/lp_data/HighsModelUtils.cpp
+++ b/highs/lp_data/HighsModelUtils.cpp
@@ -1531,6 +1531,8 @@ std::string utilPresolveRuleTypeToString(const HighsInt rule_type) {
     return "Col stuffing";
   } else if (rule_type == kPresolveRuleInitialSweep) {
     return "Initial sweep";
+  } else if (rule_type == kPresolveRuleFourierMotzkin) {
+    return "Fourier-Motzkin";
   }
   assert(1 == 0);
   return "????";
diff --git a/highs/lp_data/HighsOptions.h b/highs/lp_data/HighsOptions.h
index cfcf4463b8..f28d185e0e 100644
--- a/highs/lp_data/HighsOptions.h
+++ b/highs/lp_data/HighsOptions.h
@@ -461,6 +461,7 @@ struct HighsOptionsStruct {
   HighsInt presolve_substitution_maxfillin;
   HighsInt presolve_rule_off;
   HighsInt presolve_rule_test;
+  HighsInt presolve_fm_level;
   bool presolve_rule_logging;
   bool presolve_remove_slacks;
   bool no_unnecessary_rebuild_refactor;
@@ -636,6 +637,7 @@ struct HighsOptionsStruct {
         presolve_substitution_maxfillin(0),
         presolve_rule_off(0),
         presolve_rule_test(0),
+        presolve_fm_level(0),
         presolve_rule_logging(false),
         presolve_remove_slacks(false),
         no_unnecessary_rebuild_refactor(false),
@@ -1684,6 +1686,11 @@ class HighsOptions : public HighsOptionsStruct {
         &presolve_rule_test, 0, 0, kPresolveRuleMax);
     records.push_back(record_int);
 
+    record_int =
+        new OptionRecordInt("presolve_fm_level", "Fourier-Motzkin level",
+                            advanced, &presolve_fm_level, 0, 0, 1);
+    records.push_back(record_int);
+
     record_bool = new OptionRecordBool(
         "presolve_rule_logging", "Log effectiveness of presolve rules for LP",
         advanced, &presolve_rule_logging, false);
diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp
index be9efd1346..74380ad043 100644
--- a/highs/presolve/HPresolve.cpp
+++ b/highs/presolve/HPresolve.cpp
@@ -502,6 +502,7 @@ void HPresolve::chooseRules() {
     presolve_light_rule_off[kPresolveRuleEnumeration] = true;
     presolve_light_rule_off[kPresolveRuleDualFixing] = true;
     presolve_light_rule_off[kPresolveRuleColStuffing] = true;
+    presolve_light_rule_off[kPresolveRuleFourierMotzkin] = true;
   }
 
   if (!silent && options->log_dev_level) {
@@ -997,6 +998,8 @@ void HPresolve::shrinkProblem(HighsPostsolveStack& postsolve_stack) {
       }
     }
   }
+  if (model->fme_obj_col_ >= 0)
+    model->fme_obj_col_ = newColIndex[model->fme_obj_col_];
   colDeleted.assign(model->num_col_, false);
   model->col_cost_.resize(model->num_col_);
   model->col_lower_.resize(model->num_col_);
@@ -2409,6 +2412,7 @@ void HPresolve::markColDeleted(HighsInt col) {
     colDeleted[col] = true;
   }
   ++numDeletedCols;
+  if (col == model->fme_obj_col_) model->fme_obj_col_ = -1;
 }
 
 HPresolve::Result HPresolve::changeColUpper(HighsInt col, double newUpper) {
@@ -6331,6 +6335,8 @@ HPresolve::Result HPresolve::initialSweep(
   model->a_matrix_.start_.resize(num_col + 1);
   model->a_matrix_.index_.resize(nnz);
   model->a_matrix_.value_.resize(nnz);
+  if (model->fme_obj_col_ >= 0)
+    model->fme_obj_col_ = newColIndex[model->fme_obj_col_];
   postsolve_stack.compressColIndexMap(newColIndex);
   HPRESOLVE_CHECKED_CALL(checkLimits(postsolve_stack));
 
@@ -6717,6 +6723,7 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
         mipsolver != nullptr || !options->lp_presolve_requires_basis_postsolve;
 #endif
     bool tryProbing = mipsolver != nullptr;
+    bool tryFourierMotzkin = mipsolver != nullptr;
     HighsInt numCliquesBeforeProbing = -1;
     bool domcolAfterProbingCalled = false;
     bool dependentEquationsCalled = mipsolver != nullptr;
@@ -6753,6 +6760,13 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
             applyConflictGraphSubstitutions(postsolve_stack, numDelCol));
       }
 
+      HighsInt numColsEliminatedFourierMotzkin = 0;
+      if (tryFourierMotzkin && this->allow_rule_[kPresolveRuleFourierMotzkin]) {
+        HPRESOLVE_CHECKED_CALL(
+            fourierMotzkin(postsolve_stack, numColsEliminatedFourierMotzkin));
+        tryFourierMotzkin = false;
+      }
+
       if (reducedToEmpty()) break;
 
       if (this->allow_rule_[kPresolveRuleAggregator]) {
@@ -7772,6 +7786,864 @@ HPresolve::Result HPresolve::aggregator(HighsPostsolveStack& postsolve_stack) {
   return Result::kOk;
 }
 
+HPresolve::Result HPresolve::fourierMotzkin(
+    HighsPostsolveStack& postsolve_stack, HighsInt& numColsEliminated) {
+  assert(this->allow_rule_[kPresolveRuleFourierMotzkin]);
+  const bool logging_on = analysis_.logging_on_;
+  if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleFourierMotzkin);
+
+  using FmeRow = HighsPostsolveStack::FmeRowData<HighsTripletTreeSlicePreOrder>;
+  using FmeAncestryEntry = HighsPostsolveStack::FmeAncestryEntry;
+  using FmeBlockStep = HighsPostsolveStack::FmeBlockStep;
+
+  // max. absolute coefficient
+  const double maxCoef = 1e3;
+
+  // max. number of consecutive failures (while trying to build the heap)
+  const HighsInt maxNumFails = 100;
+  // max. size of the heap
+  const HighsInt maxHeapSize = 10000;
+
+  // sentinel row indices for variable bounds and objective row
+  const HighsInt kUpperBoundRow = -2;
+  const HighsInt kLowerBoundRow = -3;
+  const HighsInt kObjectiveRow = -4;
+
+  // structs
+  struct Heap {
+    struct Entry {
+      HighsInt col;
+      int64_t neRed;
+      int64_t mrRed;
+    };
+
+    std::vector<Entry> entries;
+    std::vector<HighsInt> pos;
+
+    bool empty() const { return entries.empty(); }
+    HighsInt size() const { return static_cast<HighsInt>(entries.size()); }
+    HighsInt top() const { return entries[0].col; }
+    bool contains(HighsInt col) const { return pos[col] != -1; }
+
+    void reset(HighsInt numCol, HighsInt reserveSize) {
+      entries.clear();
+      entries.reserve(reserveSize);
+      pos.assign(numCol, -1);
+    }
+
+    void push(HighsInt col, int64_t neRed, int64_t mrRed) {
+      pos[col] = size();
+      entries.push_back({col, neRed, mrRed});
+    }
+
+    void insert(HighsInt col, int64_t neRed, int64_t mrRed) {
+      push(col, neRed, mrRed);
+      siftUp(pos[col]);
+    }
+
+    void remove(HighsInt col) {
+      HighsInt p = pos[col];
+      if (p == -1) return;
+      swap(p, size() - 1);
+      pos[col] = -1;
+      entries.pop_back();
+      siftUp(p);
+      siftDown(p);
+    }
+
+    void update(HighsInt col, int64_t neRed, int64_t mrRed) {
+      HighsInt p = pos[col];
+      if (p == -1) return;
+      entries[p].neRed = neRed;
+      entries[p].mrRed = mrRed;
+      siftUp(p);
+      siftDown(p);
+    }
+
+    void heapify() {
+      for (HighsInt i = size() / 2 - 1; i >= 0; --i) siftDown(i);
+    }
+
+   private:
+    bool better(HighsInt i, HighsInt j) const {
+      if (entries[i].neRed != entries[j].neRed)
+        return entries[i].neRed > entries[j].neRed;
+      return entries[i].mrRed > entries[j].mrRed;
+    }
+
+    void swap(HighsInt i, HighsInt j) {
+      if (i == j) return;
+      std::swap(entries[i], entries[j]);
+      pos[entries[i].col] = i;
+      pos[entries[j].col] = j;
+    }
+
+    void siftUp(HighsInt i) {
+      if (i >= size()) return;
+      while (i > 0) {
+        HighsInt parent = (i - 1) / 2;
+        if (!better(i, parent)) break;
+        swap(i, parent);
+        i = parent;
+      }
+    }
+
+    void siftDown(HighsInt i) {
+      HighsInt n = size();
+      if (i >= n) return;
+      while (true) {
+        HighsInt best = i;
+        HighsInt left = 2 * i + 1;
+        HighsInt right = 2 * i + 2;
+        if (left < n && better(left, best)) best = left;
+        if (right < n && better(right, best)) best = right;
+        if (best == i) break;
+        swap(i, best);
+        i = best;
+      }
+    }
+  };
+
+  struct newRowEntry {
+    HighsInt col;
+    HighsCDouble val;
+  };
+
+  struct newRow {
+    std::vector<newRowEntry> entries;
+    double lower;
+    double upper;
+    HighsInt plusIndex;
+    HighsInt minusIndex;
+    double plusScale;
+    double minusScale;
+  };
+
+  struct NewRowOrigin {
+    HighsInt plusRow;
+    HighsInt minusRow;
+    double plusScale;
+    double minusScale;
+  };
+
+  auto finalise = [&]() {
+    analysis_.logging_on_ = logging_on;
+    if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleFourierMotzkin);
+    return checkLimits(postsolve_stack);
+  };
+
+  auto acceptCoef = [&](double val) {
+    double absval = std::abs(val);
+    return absval == 0.0 || (absval >= 1.0 / maxCoef && absval <= maxCoef);
+  };
+
+  auto isCandidate = [&](HighsInt col) {
+    if (colDeleted[col]) return false;
+    if (colsize[col] == 0) return false;
+    if (col == model->fme_obj_col_) return false;
+    if (model->integrality_[col] != HighsVarType::kContinuous) return false;
+    if (!acceptCoef(model->col_cost_[col])) return false;
+    if (options->presolve_fm_level < 1 && model->col_cost_[col] != 0.0)
+      return false;
+    for (const auto& nz : getColumnVector(col))
+      if (isEquation(nz.index()) || !acceptCoef(nz.value())) return false;
+    return true;
+  };
+
+  auto computeCandidates = [&](std::vector<HighsInt>& candidates) {
+    candidates.clear();
+    for (HighsInt col = 0; col < model->num_col_; col++)
+      if (isCandidate(col)) candidates.push_back(col);
+    pdqsort(candidates.begin(), candidates.end(),
+            [&](HighsInt a, HighsInt b) { return colsize[a] < colsize[b]; });
+    return !candidates.empty();
+  };
+
+  auto checkRows = [&](HighsInt col, const std::vector<HighsInt>& objRowCols,
+                       std::vector<HighsInt>& iPlus,
+                       std::vector<HighsInt>& iMinus, int64_t& nePlus,
+                       int64_t& neMinus) {
+    nePlus = 0;
+    neMinus = 0;
+    iPlus.clear();
+    iMinus.clear();
+    for (const auto& nz : getColumnVector(col)) {
+      HighsInt row = nz.index();
+      if (rowDeleted[row]) continue;
+
+      if (isRanged(row)) {
+        iPlus.push_back(row);
+        nePlus += rowsize[row];
+        iMinus.push_back(row);
+        neMinus += rowsize[row];
+      } else {
+        HighsInt direction;
+        if (model->row_lower_[row] == -kHighsInf &&
+            model->row_upper_[row] != kHighsInf)
+          direction = 1;
+        else
+          direction = -1;
+
+        if (direction * nz.value() > 0) {
+          iPlus.push_back(row);
+          nePlus += rowsize[row];
+        } else {
+          iMinus.push_back(row);
+          neMinus += rowsize[row];
+        }
+      }
+    }
+
+    // include finite variable bounds as singleton rows
+    if (model->col_upper_[col] != kHighsInf) {
+      iPlus.push_back(kUpperBoundRow);
+      nePlus += 1;
+    }
+    if (model->col_lower_[col] != -kHighsInf) {
+      iMinus.push_back(kLowerBoundRow);
+      neMinus += 1;
+    }
+
+    // simulate the objective constraint row for candidates with nonzero
+    // cost when the reformulation has not yet been performed
+    if (!objRowCols.empty() && model->col_cost_[col] != 0.0) {
+      int64_t objRowSize = static_cast<int64_t>(objRowCols.size());
+      if (model->col_cost_[col] > 0.0) {
+        iPlus.push_back(kObjectiveRow);
+        nePlus += objRowSize;
+      } else {
+        iMinus.push_back(kObjectiveRow);
+        neMinus += objRowSize;
+      }
+    }
+  };
+
+  auto collectAffectedCols = [&](HighsInt col, const std::vector<HighsInt>& set,
+                                 const std::vector<HighsInt>& objRowCols,
+                                 std::vector<HighsInt>& mark,
+                                 std::vector<HighsInt>& otherMark,
+                                 std::vector<HighsInt>& affectedCols) {
+    for (HighsInt row : set) {
+      if (row == kObjectiveRow) {
+        for (HighsInt k : objRowCols) {
+          if (k == col) continue;
+          if (mark[k] == 0 && otherMark[k] == 0) affectedCols.push_back(k);
+          mark[k]++;
+        }
+      } else {
+        if (row < 0) continue;
+        for (const auto& nz : getRowVector(row)) {
+          HighsInt k = nz.index();
+          if (k == col) continue;
+          if (mark[k] == 0 && otherMark[k] == 0) affectedCols.push_back(k);
+          mark[k]++;
+        }
+      }
+    }
+  };
+
+  auto checkNonZeros = [&](HighsInt col,
+                           const std::vector<HighsInt>& objRowCols,
+                           std::vector<HighsInt>& iPlus,
+                           std::vector<HighsInt>& iMinus,
+                           std::vector<HighsInt>& pPlus,
+                           std::vector<HighsInt>& pMinus,
+                           std::vector<HighsInt>& affectedCols, int64_t& neRed,
+                           int64_t& mrRed) {
+    // initialise
+    neRed = 0;
+    mrRed = 0;
+
+    // check rows
+    int64_t nePlus;
+    int64_t neMinus;
+    checkRows(col, objRowCols, iPlus, iMinus, nePlus, neMinus);
+
+    if (iPlus.size() == 0 || iMinus.size() == 0) {
+      // other presolve reductions may handle this case (e.g., implied free
+      // column substitution)
+      iPlus.clear();
+      iMinus.clear();
+      return false;
+    }
+
+    // take into account other variables present in the rows
+    collectAffectedCols(col, iPlus, objRowCols, pPlus, pMinus, affectedCols);
+    collectAffectedCols(col, iMinus, objRowCols, pMinus, pPlus, affectedCols);
+
+    // compute correction term
+    int64_t correction = 0;
+    for (HighsInt k : affectedCols) {
+      correction += static_cast<int64_t>(pPlus[k]) * pMinus[k];
+      pPlus[k] = 0;
+      pMinus[k] = 0;
+    }
+
+    int64_t mPlus = static_cast<int64_t>(iPlus.size());
+    int64_t mMinus = static_cast<int64_t>(iMinus.size());
+    int64_t neOld = nePlus + neMinus;
+    // note that we subtract the entries for column 'col' since these are
+    // eliminated
+    int64_t neNew =
+        mPlus * (neMinus - mMinus) + mMinus * (nePlus - mPlus) - correction;
+    neRed = neOld - neNew;
+    mrRed = mPlus + mMinus - mPlus * mMinus;
+    return true;
+  };
+
+  auto checkNewRow = [&](const newRow& nr, bool& isRedundant) {
+    HighsCDouble impliedLower = 0;
+    HighsCDouble impliedUpper = 0;
+    bool lowerFinite = true;
+    bool upperFinite = true;
+    isRedundant = false;
+    for (const auto& e : nr.entries) {
+      double lb = model->col_lower_[e.col];
+      double ub = model->col_upper_[e.col];
+      if (e.val > 0) {
+        lowerFinite = lowerFinite && lb != -kHighsInf;
+        if (lowerFinite) impliedLower += e.val * lb;
+        upperFinite = upperFinite && ub != kHighsInf;
+        if (upperFinite) impliedUpper += e.val * ub;
+      } else {
+        lowerFinite = lowerFinite && ub != kHighsInf;
+        if (lowerFinite) impliedLower += e.val * ub;
+        upperFinite = upperFinite && lb != -kHighsInf;
+        if (upperFinite) impliedUpper += e.val * lb;
+      }
+      if (!lowerFinite && !upperFinite) return Result::kOk;
+    }
+
+    double lower = lowerFinite ? static_cast<double>(impliedLower) : -kHighsInf;
+    double upper = upperFinite ? static_cast<double>(impliedUpper) : kHighsInf;
+
+    // check for infeasibility
+    if (lower > nr.upper + primal_feastol || upper < nr.lower - primal_feastol)
+      return Result::kPrimalInfeasible;
+
+    // check for redundancy
+    isRedundant = lower >= nr.lower - primal_feastol &&
+                  upper <= nr.upper + primal_feastol;
+
+    return Result::kOk;
+  };
+
+  auto getRowData = [&](HighsInt row, HighsInt col, HighsInt multiplier,
+                        double& absCoef, HighsInt& direction, double& bound) {
+    if (row < 0) {
+      // artificial lower / upper bound row
+      direction = 1;
+      absCoef = 1.0;
+      bound = multiplier > 0 ? model->col_upper_[col] : -model->col_lower_[col];
+    } else {
+      HighsInt pPos = findNonzero(row, col);
+      assert(pPos != -1);
+      direction = multiplier * Avalue[pPos] > 0 ? HighsInt{1} : HighsInt{-1};
+      absCoef = std::abs(Avalue[pPos]);
+      bound = direction > 0 ? model->row_upper_[row] : -model->row_lower_[row];
+    }
+  };
+
+  auto collectRowEntries = [&](HighsInt row, HighsInt col, double scale,
+                               std::vector<newRowEntry>& newRowEntries,
+                               std::vector<HighsInt>& newRowMark) {
+    if (row < 0) return;
+    for (const auto& nz : getRowVector(row)) {
+      if (nz.index() == col) continue;
+      double val = scale * nz.value();
+      if (newRowMark[nz.index()] == -1) {
+        newRowMark[nz.index()] = static_cast<HighsInt>(newRowEntries.size());
+        newRowEntries.push_back({nz.index(), val});
+      } else {
+        newRowEntries[newRowMark[nz.index()]].val += val;
+      }
+    }
+  };
+
+  auto isReduction = [](int64_t neRed, int64_t mrRed) {
+    return neRed > 0 || (neRed == 0 && mrRed > 0);
+  };
+
+  auto insertOriginals =
+      [&](std::set<HighsInt>& rows,
+          const std::unordered_map<HighsInt, std::set<HighsInt>>& originals,
+          HighsInt row, HighsInt col) {
+        if (row == kUpperBoundRow)
+          rows.insert(-(2 * col + 1));
+        else if (row == kLowerBoundRow)
+          rows.insert(-(2 * col + 2));
+        else {
+          auto it = originals.find(row);
+          if (it != originals.end())
+            rows.insert(it->second.begin(), it->second.end());
+          else
+            rows.insert(row);
+        }
+      };
+
+  auto mergeOriginals =
+      [&](std::set<HighsInt>& rows,
+          const std::unordered_map<HighsInt, std::set<HighsInt>>& originals,
+          HighsInt plusRow, HighsInt minusRow, HighsInt col) {
+        rows.clear();
+        insertOriginals(rows, originals, plusRow, col);
+        insertOriginals(rows, originals, minusRow, col);
+      };
+
+  auto cernikovRedundant =
+      [&](std::set<HighsInt>& rows,
+          const std::unordered_map<HighsInt, std::set<HighsInt>>& originals,
+          HighsInt plusRow, HighsInt minusRow, HighsInt col,
+          HighsInt numColsElim) {
+        mergeOriginals(rows, originals, plusRow, minusRow, col);
+        return static_cast<HighsInt>(rows.size()) > numColsElim + 2;
+      };
+
+  // reformulate objective as a constraint: min c^T x + offset becomes
+  // min z with c^T x - z <= -offset. this allows FME to eliminate
+  // continuous columns with nonzero cost.
+  auto reformulateObjective = [&]() {
+    if (model->fme_obj_col_ != -1) {
+      assert(!colDeleted[model->fme_obj_col_]);
+      return;
+    }
+
+    HighsInt zCol = model->num_col_;
+    model->num_col_++;
+    model->a_matrix_.num_col_++;
+
+    // extend model vectors
+    model->col_cost_.push_back(1.0);
+    model->col_lower_.push_back(-kHighsInf);
+    model->col_upper_.push_back(kHighsInf);
+    model->integrality_.push_back(HighsVarType::kContinuous);
+    model->a_matrix_.start_.push_back(model->a_matrix_.start_.back());
+    if (model->col_names_.size() > 0) model->col_names_.push_back("fme_obj_z");
+
+    // extend presolve vectors
+    colhead.push_back(-1);
+    colsize.push_back(0);
+    colDeleted.push_back(0);
+    implColLower.push_back(-kHighsInf);
+    implColUpper.push_back(kHighsInf);
+    colLowerSource.push_back(-1);
+    colUpperSource.push_back(-1);
+    implRowDualSourceByCol.push_back({});
+    changedColFlag.push_back(1);
+    numProbes.push_back(0);
+
+    // update implied bound structures (pointers may be invalidated by
+    // reallocation of column vectors above)
+    impliedRowBounds.setBoundArrays(
+        model->col_lower_.data(), model->col_upper_.data(), implColLower.data(),
+        implColUpper.data(), colLowerSource.data(), colUpperSource.data());
+    impliedDualRowBounds.setNumSums(model->num_col_);
+
+    // register in postsolve stack
+    postsolve_stack.appendColToModel();
+
+    // build the objective constraint row: c^T x - z <= -offset
+    double offset = model->offset_;
+    std::vector<HighsInt> objIndices;
+    std::vector<double> objValues;
+    for (HighsInt j = 0; j < zCol; ++j) {
+      if (!colDeleted[j] && model->col_cost_[j] != 0.0) {
+        objIndices.push_back(j);
+        objValues.push_back(model->col_cost_[j]);
+      }
+    }
+    objIndices.push_back(zCol);
+    objValues.push_back(-1.0);
+
+    // zero out original costs and offset
+    for (HighsInt j = 0; j < zCol; ++j) model->col_cost_[j] = 0.0;
+    model->offset_ = 0.0;
+
+    // add the constraint row to the matrix
+    addToMatrix(postsolve_stack, -kHighsInf, -offset, objIndices, objValues);
+
+    // register reduction so getReducedPrimalSolution can compute z
+    std::vector<HighsPostsolveStack::Nonzero> costEntries;
+    for (size_t k = 0; k < objIndices.size(); k++)
+      costEntries.emplace_back(objIndices[k], objValues[k]);
+    postsolve_stack.fourierMotzkinObjCol(zCol, offset, costEntries);
+
+    model->fme_obj_col_ = zCol;
+
+    shrinkProblem(postsolve_stack);
+  };
+
+  auto collectCandidatesAndBuildHeap =
+      [&](std::vector<HighsInt>& candidates, Heap& heap,
+          std::vector<HighsInt>& iPlus, std::vector<HighsInt>& iMinus,
+          std::vector<HighsInt>& pPlus, std::vector<HighsInt>& pMinus,
+          std::vector<HighsInt>& affectedCols,
+          const std::vector<HighsInt>& objRowCols) {
+        // compute candidates
+        if (!computeCandidates(candidates)) return false;
+        // set up data structures for heap
+        heap.reset(model->num_col_, static_cast<HighsInt>(candidates.size()));
+        pPlus.assign(model->num_col_, 0);
+        pMinus.assign(model->num_col_, 0);
+        iPlus.reserve(model->num_row_);
+        iMinus.reserve(model->num_row_);
+        affectedCols.reserve(model->num_col_);
+        // inspect candidates (with limits)
+        HighsInt numFails = 0;
+        for (HighsInt col : candidates) {
+          int64_t neRed;
+          int64_t mrRed;
+          bool elimCandidate =
+              checkNonZeros(col, objRowCols, iPlus, iMinus, pPlus, pMinus,
+                            affectedCols, neRed, mrRed);
+          affectedCols.clear();
+          if (!elimCandidate || !isReduction(neRed, mrRed)) {
+            // count number of failures
+            if (++numFails > maxNumFails) break;
+            continue;
+          }
+          // add to heap
+          numFails = 0;
+          heap.push(col, neRed, mrRed);
+          if (heap.size() >= maxHeapSize) break;
+        }
+        if (heap.empty()) return false;
+        // heapify
+        heap.heapify();
+        return true;
+      };
+
+  // find index of a row within a list
+  auto findRowIndex = [](HighsInt row,
+                         const std::vector<FmeRow>& rows) -> HighsInt {
+    for (HighsInt i = 0; i < static_cast<HighsInt>(rows.size()); ++i)
+      if (rows[i].row == row) return i;
+    return -1;
+  };
+
+  auto collectRows = [&](const std::vector<HighsInt>& rows) {
+    std::vector<FmeRow> result;
+    for (HighsInt r : rows) {
+      if (r < 0) continue;
+      result.push_back(
+          {r, model->row_lower_[r], model->row_upper_[r], getRowVector(r)});
+    }
+    return result;
+  };
+
+  auto inheritAncestry =
+      [&](std::unordered_map<HighsInt, std::vector<FmeAncestryEntry>>&
+              rowAncestry,
+          HighsInt newModelRow, HighsInt parentRow, HighsInt parentRowIndex,
+          HighsInt stepIndex, double scale, bool isMinus) {
+        if (parentRow < 0) return;
+        auto it = rowAncestry.find(parentRow);
+        if (it != rowAncestry.end()) {
+          for (const auto& a : it->second)
+            rowAncestry[newModelRow].push_back(
+                {a.step, a.parentRowIndex, a.scale * scale, a.isMinus});
+        }
+        if (parentRowIndex >= 0)
+          rowAncestry[newModelRow].push_back(
+              {stepIndex, parentRowIndex, scale, isMinus});
+      };
+
+  auto printLog = [&](HighsInt colsRemoved, HighsInt rowsRemoved,
+                      HighsInt rowsAdded) {
+    highsLogDev(options->log_options, HighsLogType::kInfo,
+                "Fourier-Motzkin (%s objective reformulation) added "
+                "%" HIGHSINT_FORMAT " rows and eliminated %" HIGHSINT_FORMAT
+                " rows and %" HIGHSINT_FORMAT " columns\n",
+                options->presolve_fm_level >= 1 ? "with" : "without", rowsAdded,
+                rowsRemoved, colsRemoved);
+  };
+
+  // workspace vectors
+  std::vector<HighsInt> candidates;
+  std::vector<HighsInt> iPlus;
+  std::vector<HighsInt> iMinus;
+  std::vector<HighsInt> pPlus;
+  std::vector<HighsInt> pMinus;
+  std::vector<HighsInt> affectedCols;
+
+  // indexed max-heap
+  Heap heap;
+
+  // precompute the objective row: columns with nonzero cost
+  // used to simulate the objective constraint in checkRows before
+  // reformulation actually happens
+  std::vector<HighsInt> objRowCols;
+  if (model->fme_obj_col_ == -1 && options->presolve_fm_level >= 1) {
+    for (HighsInt j = 0; j < model->num_col_; ++j) {
+      if (!colDeleted[j] && model->col_cost_[j] != 0.0) objRowCols.push_back(j);
+    }
+  }
+
+  // compute candidates and build initial heap
+  if (!collectCandidatesAndBuildHeap(candidates, heap, iPlus, iMinus, pPlus,
+                                     pMinus, affectedCols, objRowCols))
+    return finalise();
+
+  // vectors for computing new row entries
+  std::vector<newRowEntry> newRowEntries;
+  std::vector<HighsInt> newRowMark(model->num_col_, -1);
+
+  // vector for storing new rows
+  std::vector<newRow> newRows;
+
+  // workspace for filtering new rows
+  std::vector<double> rowLower;
+  std::vector<double> rowUpper;
+  std::vector<std::vector<HighsInt>> rowIndices;
+  std::vector<std::vector<double>> rowValues;
+  std::vector<NewRowOrigin> newRowOrigins;
+
+  // vector for saving affected candidates
+  std::vector<HighsInt> saveAffectedCols;
+
+  // counters for numbers of eliminations
+  numColsEliminated = 0;
+  HighsInt numColsEliminatedBlock = 0;
+  HighsInt numRowsEliminated = 0;
+  HighsInt numRowsAdded = 0;
+
+  // FM block data for postsolve
+  std::vector<FmeBlockStep> blockSteps;
+
+  // surviving row to its ancestry (which parent rows it descends from)
+  std::unordered_map<HighsInt, std::vector<FmeAncestryEntry>> rowAncestry;
+
+  // distinct original parent rows for each derived row (Cernikov check)
+  std::unordered_map<HighsInt, std::set<HighsInt>> rowOriginals;
+  std::set<HighsInt> mergedOriginals;
+
+  // main loop: eliminate variables from heap
+  while (!heap.empty()) {
+    HighsInt col = heap.top();
+    heap.remove(col);
+
+    // if this candidate has nonzero cost and objective has not yet been
+    // reformulated, perform the reformulation now and rebuild the heap
+    if (model->fme_obj_col_ == -1 && model->col_cost_[col] != 0.0) {
+      // finalise any in-progress FM block before reformulating, since
+      // reformulateObjective pushes other reductions onto the data stack
+      if (!blockSteps.empty()) {
+        postsolve_stack.fourierMotzkinBlockFinalise(blockSteps, rowAncestry);
+        printLog(numColsEliminatedBlock, numRowsEliminated, numRowsAdded);
+        blockSteps.clear();
+        rowAncestry.clear();
+        rowOriginals.clear();
+        numColsEliminatedBlock = 0;
+        numRowsEliminated = 0;
+        numRowsAdded = 0;
+      }
+      // reformulate objective
+      reformulateObjective();
+      // clear vector for objective and resize marker
+      objRowCols.clear();
+      newRowMark.resize(model->num_col_, -1);
+      // re-compute candidates and re-build heap
+      if (!collectCandidatesAndBuildHeap(candidates, heap, iPlus, iMinus, pPlus,
+                                         pMinus, affectedCols, objRowCols))
+        return finalise();
+      continue;
+    }
+
+    // compute affected columns
+    int64_t neRed;
+    int64_t mrRed;
+    bool elimCandidate = checkNonZeros(col, objRowCols, iPlus, iMinus, pPlus,
+                                       pMinus, affectedCols, neRed, mrRed);
+
+    // heap data should be up-to-date
+    assert(elimCandidate && isReduction(neRed, mrRed));
+
+    HighsInt stepIdx = static_cast<HighsInt>(blockSteps.size());
+
+    // perform elimination: generate new rows
+    newRows.clear();
+    for (HighsInt pRow : iPlus) {
+      double pCoefAbs;
+      double pBound;
+      HighsInt pDirection;
+      getRowData(pRow, col, HighsInt{1}, pCoefAbs, pDirection, pBound);
+
+      for (HighsInt mRow : iMinus) {
+        double mCoefAbs;
+        double mBound;
+        HighsInt mDirection;
+        getRowData(mRow, col, HighsInt{-1}, mCoefAbs, mDirection, mBound);
+
+        // scale factor to preserve violation tolerances (see section 4.3):
+        double s = (pCoefAbs * mCoefAbs) / (pCoefAbs + mCoefAbs);
+        double pScale = s / pCoefAbs;
+        double mScale = s / mCoefAbs;
+
+        // collect row entries
+        collectRowEntries(pRow, col, pDirection * pScale, newRowEntries,
+                          newRowMark);
+        collectRowEntries(mRow, col, mDirection * mScale, newRowEntries,
+                          newRowMark);
+
+        // reset marker before removing near-zeros
+        for (const auto& e : newRowEntries) newRowMark[e.col] = -1;
+
+        // remove near-zero entries
+        newRowEntries.erase(
+            std::remove_if(newRowEntries.begin(), newRowEntries.end(),
+                           [&](const newRowEntry& e) {
+                             return abs(e.val) <= options->small_matrix_value;
+                           }),
+            newRowEntries.end());
+
+        // store new row
+        double new_upper =
+            static_cast<double>(static_cast<HighsCDouble>(pScale) * pBound +
+                                static_cast<HighsCDouble>(mScale) * mBound);
+        newRows.push_back({newRowEntries, -kHighsInf, new_upper, pRow, mRow,
+                           pDirection * pScale, mDirection * mScale});
+
+        // clear vector
+        newRowEntries.clear();
+      }
+    }
+
+    // add new rows, filtering out redundant ones
+    rowLower.clear();
+    rowUpper.clear();
+    rowIndices.clear();
+    rowValues.clear();
+    newRowOrigins.clear();
+
+    for (const auto& nr : newRows) {
+      bool redundant = false;
+      HPRESOLVE_CHECKED_CALL(checkNewRow(nr, redundant));
+      if (redundant) continue;
+
+      // Cernikov redundancy check
+      if (cernikovRedundant(mergedOriginals, rowOriginals, nr.plusIndex,
+                            nr.minusIndex, col, numColsEliminated))
+        continue;
+
+      std::vector<HighsInt> indices;
+      std::vector<double> values;
+      indices.reserve(nr.entries.size());
+      values.reserve(nr.entries.size());
+      for (const auto& e : nr.entries) {
+        indices.push_back(e.col);
+        values.push_back(static_cast<double>(e.val));
+      }
+      rowLower.push_back(nr.lower);
+      rowUpper.push_back(nr.upper);
+      rowIndices.push_back(std::move(indices));
+      rowValues.push_back(std::move(values));
+      newRowOrigins.push_back(
+          {nr.plusIndex, nr.minusIndex, nr.plusScale, nr.minusScale});
+    }
+
+    // serialize row data for postsolve before addToMatrix invalidates slices
+    std::vector<FmeRow> plusRows = collectRows(iPlus);
+    std::vector<FmeRow> minusRows = collectRows(iMinus);
+
+    // push row data for this elimination step onto the postsolve stack
+    postsolve_stack.fourierMotzkinBlockPushStep(col, plusRows, minusRows);
+
+    // save block metadata
+    assert(model->col_cost_[col] == 0.0);
+    blockSteps.push_back({col,
+                          model->col_lower_[col],
+                          model->col_upper_[col],
+                          static_cast<HighsInt>(plusRows.size()),
+                          static_cast<HighsInt>(minusRows.size()),
+                          {}});
+
+    // add new rows to matrix
+    HighsInt firstNewRow = model->num_row_;
+    if (!addToMatrix(postsolve_stack, rowLower, rowUpper, rowIndices,
+                     rowValues))
+      return finalise();
+    numRowsAdded += static_cast<HighsInt>(rowIndices.size());
+
+    // build FmeNewRow data and ancestry for this step
+    auto& stepNewRows = blockSteps.back().newRows;
+    stepNewRows.reserve(newRowOrigins.size());
+    for (HighsInt k = 0; k < static_cast<HighsInt>(newRowOrigins.size()); ++k) {
+      HighsInt newModelRow = firstNewRow + k;
+      const auto& origin = newRowOrigins[k];
+      HighsInt pIdx = findRowIndex(origin.plusRow, plusRows);
+      HighsInt mIdx = findRowIndex(origin.minusRow, minusRows);
+      inheritAncestry(rowAncestry, newModelRow, origin.plusRow, pIdx, stepIdx,
+                      origin.plusScale, false);
+      inheritAncestry(rowAncestry, newModelRow, origin.minusRow, mIdx, stepIdx,
+                      origin.minusScale, true);
+      mergeOriginals(mergedOriginals, rowOriginals, origin.plusRow,
+                     origin.minusRow, col);
+      rowOriginals[newModelRow] = mergedOriginals;
+      stepNewRows.push_back({newModelRow, pIdx, mIdx});
+    }
+
+    // mark column as deleted
+    markColDeleted(col);
+    ++numColsEliminatedBlock;
+    ++numColsEliminated;
+
+    // remove old rows containing col (skip bound rows)
+    for (HighsInt rp : iPlus) {
+      if (rp < 0) continue;
+      rowAncestry.erase(rp);
+      rowOriginals.erase(rp);
+      removeRow(rp);
+      ++numRowsEliminated;
+    }
+    for (HighsInt rm : iMinus) {
+      if (rm < 0) continue;
+      rowAncestry.erase(rm);
+      rowOriginals.erase(rm);
+      if (rowDeleted[rm]) continue;
+      removeRow(rm);
+      ++numRowsEliminated;
+    }
+
+    // update affected candidates in the heap
+    saveAffectedCols.swap(affectedCols);
+    for (HighsInt k : saveAffectedCols) {
+      // check if variable is a candidate
+      bool isCandidateCol = isCandidate(k);
+      // skip variable if it is not on the heap and no candidate
+      if (!heap.contains(k) && !isCandidateCol) continue;
+      // check column non-zeros
+      int64_t ne, mr;
+      bool elimCandidate =
+          isCandidateCol && checkNonZeros(k, objRowCols, iPlus, iMinus, pPlus,
+                                          pMinus, affectedCols, ne, mr);
+      affectedCols.clear();
+      if (!elimCandidate || !isReduction(ne, mr)) {
+        // no candidate or not beneficial -> remove from heap
+        heap.remove(k);
+      } else if (!heap.contains(k)) {
+        // new candidate -> insert into heap
+        heap.insert(k, ne, mr);
+      } else {
+        // update heap
+        heap.update(k, ne, mr);
+      }
+    }
+    saveAffectedCols.clear();
+
+    if (checkLimits(postsolve_stack) != Result::kOk) break;
+  }
+
+  if (numColsEliminatedBlock > 0) {
+    // finalize the FM block
+    postsolve_stack.fourierMotzkinBlockFinalise(blockSteps, rowAncestry);
+
+    // log message
+    printLog(numColsEliminatedBlock, numRowsEliminated, numRowsAdded);
+  }
+
+  return finalise();
+}
+
 void HPresolve::substitute(HighsInt substcol, HighsInt staycol, double offset,
                            double scale) {
   // substitute the column in each row where it occurs
diff --git a/highs/presolve/HPresolve.h b/highs/presolve/HPresolve.h
index 871a69a223..9e6884017f 100644
--- a/highs/presolve/HPresolve.h
+++ b/highs/presolve/HPresolve.h
@@ -505,6 +505,9 @@ class HPresolve {
 
   Result aggregator(HighsPostsolveStack& postsolve_stack);
 
+  Result fourierMotzkin(HighsPostsolveStack& postsolve_stack,
+                        HighsInt& numColsEliminated);
+
   Result removeRowSingletons(HighsPostsolveStack& postsolve_stack);
 
   Result presolveColSingletons(HighsPostsolveStack& postsolve_stack);
@@ -548,6 +551,7 @@ class HPresolve {
   Result presolveRuleTestColStuffing(HighsPostsolveStack& postsolve_stack);
   Result presolveRuleTestParallelRowsAndCols(
       HighsPostsolveStack& postsolve_stack);
+  Result presolveRuleTestFourierMotzkin(HighsPostsolveStack& postsolve_stack);
 
   // Not currently called
   static void debug(const HighsLp& lp, const HighsOptions& options);
diff --git a/highs/presolve/HPresolveTest.cpp b/highs/presolve/HPresolveTest.cpp
index ff06d91285..3f248ef37a 100644
--- a/highs/presolve/HPresolveTest.cpp
+++ b/highs/presolve/HPresolveTest.cpp
@@ -16,6 +16,8 @@ HPresolve::Result HPresolve::presolveRuleTest(
     return presolveRuleTestColStuffing(postsolve_stack);
   } else if (options->presolve_rule_test == kPresolveRuleParallelRowsAndCols) {
     return presolveRuleTestParallelRowsAndCols(postsolve_stack);
+  } else if (options->presolve_rule_test == kPresolveRuleFourierMotzkin) {
+    return presolveRuleTestFourierMotzkin(postsolve_stack);
   }
   return Result::kOk;
 }
@@ -34,7 +36,8 @@ HPresolve::Result HPresolve::presolveRuleTestColStuffing(
   highsLogUser(options->log_options, HighsLogType::kInfo,
                "HPresolve::presolveRuleTestColStuffing: Stuffing removed %d "
                "rows and %d columns\n",
-               int(numDeletedRows), int(numDeletedCols));
+               static_cast<int>(numDeletedRows),
+               static_cast<int>(numDeletedCols));
   // Possibly remove the row
   return rowPresolve(postsolve_stack, 0);
 }
@@ -45,4 +48,21 @@ HPresolve::Result HPresolve::presolveRuleTestParallelRowsAndCols(
                "HPresolve::presolveRuleTestParallelRowsAndCols\n");
   return detectParallelRowsAndCols(postsolve_stack);
 }
+HPresolve::Result HPresolve::presolveRuleTestFourierMotzkin(
+    HighsPostsolveStack& postsolve_stack) {
+  assert(options->presolve_rule_test == kPresolveRuleFourierMotzkin);
+  highsLogUser(options->log_options, HighsLogType::kInfo,
+               "HPresolve::presolveRuleTestFourierMotzkin\n");
+
+  HighsInt numColsEliminated;
+  HPresolve::Result result = fourierMotzkin(postsolve_stack, numColsEliminated);
+  if (result != Result::kOk) return result;
+
+  highsLogUser(options->log_options, HighsLogType::kInfo,
+               "HPresolve::presolveRuleTestFourierMotzkin: Removed %d "
+               "rows and %d columns\n",
+               static_cast<int>(numDeletedRows),
+               static_cast<int>(numDeletedCols));
+  return result;
+}
 }  // namespace presolve
diff --git a/highs/presolve/HighsPostsolveStack.cpp b/highs/presolve/HighsPostsolveStack.cpp
index 317cd3ff4d..bc3c19600f 100644
--- a/highs/presolve/HighsPostsolveStack.cpp
+++ b/highs/presolve/HighsPostsolveStack.cpp
@@ -1345,6 +1345,386 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
   }
 }
 
+void HighsPostsolveStack::FourierMotzkinObjCol::transformToPresolvedSpace(
+    const std::vector<Nonzero>& costEntries,
+    std::vector<double>& primalSol) const {
+  double val = offset;
+  for (const Nonzero& entry : costEntries)
+    val += entry.value * primalSol[entry.index];
+  primalSol[col] = val;
+}
+
+void HighsPostsolveStack::FourierMotzkinObjCol::undo(
+    const std::vector<Nonzero>& costEntries, HighsSolution& solution) const {
+  if (!solution.dual_valid) return;
+  double zDual = solution.col_dual[col];
+  for (const Nonzero& entry : costEntries)
+    solution.col_dual[entry.index] += entry.value * zDual;
+  solution.col_dual[col] = 0.0;
+}
+
+std::vector<HighsPostsolveStack::FmeStepData>
+HighsPostsolveStack::popFourierMotzkinBlock(HighsDataStack& stack) {
+  HighsInt numSteps;
+  stack.pop(numSteps);
+
+  std::vector<FmeStepData> steps(numSteps);
+
+  // step headers
+  for (HighsInt s = numSteps - 1; s >= 0; --s) stack.pop(steps[s].header);
+
+  // new row origins
+  for (HighsInt s = numSteps - 1; s >= 0; --s) stack.pop(steps[s].newRows);
+
+  // descendants
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    HighsInt numMinus = steps[s].header.numMinus;
+    steps[s].minusDescendants.resize(numMinus);
+    for (HighsInt m = numMinus - 1; m >= 0; --m)
+      stack.pop(steps[s].minusDescendants[m]);
+    HighsInt numPlus = steps[s].header.numPlus;
+    steps[s].plusDescendants.resize(numPlus);
+    for (HighsInt p = numPlus - 1; p >= 0; --p)
+      stack.pop(steps[s].plusDescendants[p]);
+  }
+
+  // row data
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    // minus row data
+    stack.pop(steps[s].minusHeaders);
+    stack.pop(steps[s].minusCoefs);
+    HighsInt numMinus = static_cast<HighsInt>(steps[s].minusCoefs.size());
+    steps[s].minusEntries.resize(numMinus);
+    for (HighsInt r = numMinus - 1; r >= 0; --r)
+      stack.pop(steps[s].minusEntries[r]);
+
+    // plus row data
+    stack.pop(steps[s].plusHeaders);
+    stack.pop(steps[s].plusCoefs);
+    HighsInt numPlus = static_cast<HighsInt>(steps[s].plusCoefs.size());
+    steps[s].plusEntries.resize(numPlus);
+    for (HighsInt r = numPlus - 1; r >= 0; --r)
+      stack.pop(steps[s].plusEntries[r]);
+  }
+
+  return steps;
+}
+
+void HighsPostsolveStack::undoFourierMotzkinBlock(
+    const std::vector<FmeStepData>& steps, const HighsOptions& options,
+    HighsSolution& solution, HighsBasis& basis) {
+  const double tol = options.primal_feasibility_tolerance;
+  const double dual_tol = options.dual_feasibility_tolerance;
+
+  HighsInt numSteps = static_cast<HighsInt>(steps.size());
+
+  // primal postsolve (Algorithm 3): process in reverse elimination order
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    const auto& step = steps[s];
+    HighsInt col = step.header.col;
+    double lower = step.header.colLower;
+    double upper = step.header.colUpper;
+
+    auto tightenBounds = [&](const std::vector<FmeRowHeader>& headers,
+                             const std::vector<double>& coefs,
+                             const std::vector<std::vector<Nonzero>>& entries,
+                             double& lowerBound, double& upperBound) {
+      for (size_t r = 0; r < headers.size(); ++r) {
+        double aij = coefs[r];
+        HighsCDouble sum = 0.0;
+        for (const auto& nz : entries[r])
+          sum += static_cast<HighsCDouble>(nz.value) *
+                 solution.col_value[nz.index];
+        HighsInt direction = aij > 0 ? HighsInt{1} : HighsInt{-1};
+        double rhs_upper =
+            direction > 0 ? headers[r].rowUpper : headers[r].rowLower;
+        double rhs_lower =
+            direction > 0 ? headers[r].rowLower : headers[r].rowUpper;
+        if (direction * rhs_upper != kHighsInf) {
+          double bound = static_cast<double>((rhs_upper - sum) / aij);
+          upperBound = std::min(upperBound, bound);
+        }
+        if (direction * rhs_lower != -kHighsInf) {
+          double bound = static_cast<double>((rhs_lower - sum) / aij);
+          lowerBound = std::max(lowerBound, bound);
+        }
+      }
+    };
+
+    tightenBounds(step.plusHeaders, step.plusCoefs, step.plusEntries, lower,
+                  upper);
+    tightenBounds(step.minusHeaders, step.minusCoefs, step.minusEntries, lower,
+                  upper);
+
+    if (lower <= tol && upper >= -tol)
+      solution.col_value[col] = 0.0;
+    else if (lower > 0.0)
+      solution.col_value[col] = lower;
+    else
+      solution.col_value[col] = upper;
+  }
+
+  if (!solution.dual_valid) return;
+
+  // dual postsolve (Algorithm 4): process in reverse elimination order
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    const auto& step = steps[s];
+    HighsInt col = step.header.col;
+    HighsInt numPlus = step.header.numPlus;
+    HighsInt numMinus = step.header.numMinus;
+
+    // u_i = Σ_{k ∈ K^j_i} λ_k * scaleFactor
+    auto recoverDual =
+        [&](const std::vector<FmeRowHeader>& headers,
+            const std::vector<std::vector<FmeDescendant>>& descendants) {
+          for (size_t r = 0; r < headers.size(); ++r) {
+            HighsCDouble dual = 0.0;
+            for (const auto& desc : descendants[r])
+              dual += static_cast<HighsCDouble>(solution.row_dual[desc.row]) *
+                      desc.scaleFactor;
+            solution.row_dual[headers[r].row] += static_cast<double>(dual);
+          }
+        };
+    recoverDual(step.plusHeaders, step.plusDescendants);
+    recoverDual(step.minusHeaders, step.minusDescendants);
+
+    // col_dual = -Σ a_{ij} * row_dual[i] (cost is zero after reformulation)
+    HighsCDouble colDual = 0.0;
+    std::vector<HighsBool> visited(solution.row_dual.size(), false);
+    for (HighsInt r = 0; r < numPlus; ++r) {
+      HighsInt row = step.plusHeaders[r].row;
+      colDual -=
+          static_cast<HighsCDouble>(step.plusCoefs[r]) * solution.row_dual[row];
+      visited[row] = true;
+    }
+    for (HighsInt r = 0; r < numMinus; ++r) {
+      HighsInt row = step.minusHeaders[r].row;
+      if (visited[row]) continue;
+      colDual -= static_cast<HighsCDouble>(step.minusCoefs[r]) *
+                 solution.row_dual[row];
+    }
+    solution.col_dual[col] = static_cast<double>(colDual);
+  }
+
+  // basis postsolve: use dual solution to determine basis status
+  if (!basis.valid) return;
+
+  // pre-compute lower and upper slacks for each row
+  auto computeSlacks =
+      [&](HighsInt col, const std::vector<FmeRowHeader>& headers,
+          const std::vector<double>& coefs,
+          const std::vector<std::vector<Nonzero>>& entries,
+          std::vector<double>& lowerSlacks, std::vector<double>& upperSlacks) {
+        HighsInt n = static_cast<HighsInt>(headers.size());
+        lowerSlacks.resize(n);
+        upperSlacks.resize(n);
+        for (HighsInt r = 0; r < n; ++r) {
+          HighsCDouble activity =
+              static_cast<HighsCDouble>(coefs[r]) * solution.col_value[col];
+          for (const auto& nz : entries[r])
+            activity += static_cast<HighsCDouble>(nz.value) *
+                        solution.col_value[nz.index];
+          double act = static_cast<double>(activity);
+          lowerSlacks[r] = headers[r].rowLower != -kHighsInf
+                               ? act - headers[r].rowLower
+                               : kHighsInf;
+          upperSlacks[r] = headers[r].rowUpper != kHighsInf
+                               ? headers[r].rowUpper - act
+                               : kHighsInf;
+        }
+      };
+
+  // row must be basic if it has zero dual and activity strictly
+  // between bounds (complementary slackness)
+  auto rowMustBeBasic = [&](HighsInt row, double lowerSlack,
+                            double upperSlack) {
+    return std::abs(solution.row_dual[row]) <= dual_tol && lowerSlack > tol &&
+           upperSlack > tol;
+  };
+
+  // assign row as basic
+  auto assignBasicRowStatus = [&](HighsInt row, HighsInt& basicAssigned) {
+    if (basis.row_status[row] == HighsBasisStatus::kBasic ||
+        std::abs(solution.row_dual[row]) > dual_tol)
+      return false;
+    basis.row_status[row] = HighsBasisStatus::kBasic;
+    basicAssigned++;
+    return true;
+  };
+
+  // assign row as non-basic
+  auto assignNonBasicRowStatus = [&](HighsInt row, double lowerSlack,
+                                     double upperSlack) {
+    if (solution.row_dual[row] > dual_tol)
+      basis.row_status[row] = HighsBasisStatus::kLower;
+    else if (solution.row_dual[row] < -dual_tol)
+      basis.row_status[row] = HighsBasisStatus::kUpper;
+    else
+      basis.row_status[row] = upperSlack < lowerSlack
+                                  ? HighsBasisStatus::kUpper
+                                  : HighsBasisStatus::kLower;
+  };
+
+  // assign row status
+  auto assignRowStatus = [&](HighsInt row, double lowerSlack, double upperSlack,
+                             HighsInt& basicAssigned,
+                             bool forceNonBasic = false) {
+    if (forceNonBasic || !assignBasicRowStatus(row, basicAssigned))
+      assignNonBasicRowStatus(row, lowerSlack, upperSlack);
+  };
+
+  // collect a single candidate for basic assignment
+  auto collectCandidate =
+      [&](HighsInt row, bool forcedNonBasic, double lowerSlack,
+          double upperSlack, const std::vector<Nonzero>& entries,
+          HighsInt parentIndex, bool isMinus,
+          std::vector<std::tuple<HighsInt, HighsInt, bool>>& candidates) {
+        if (basis.row_status[row] == HighsBasisStatus::kBasic) return;
+        if (forcedNonBasic || std::abs(solution.row_dual[row]) > dual_tol) {
+          assignNonBasicRowStatus(row, lowerSlack, upperSlack);
+          return;
+        }
+        HighsInt nonBasicCount = 0;
+        for (const auto& nz : entries)
+          if (basis.col_status[nz.index] != HighsBasisStatus::kBasic)
+            nonBasicCount++;
+        candidates.emplace_back(nonBasicCount, parentIndex, isMinus);
+      };
+
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    const auto& step = steps[s];
+    HighsInt col = step.header.col;
+    HighsInt numPlus = step.header.numPlus;
+    HighsInt numMinus = step.header.numMinus;
+
+    // compute slacks
+    std::vector<double> plusLowerSlack;
+    std::vector<double> plusUpperSlack;
+    std::vector<double> minusLowerSlack;
+    std::vector<double> minusUpperSlack;
+    computeSlacks(col, step.plusHeaders, step.plusCoefs, step.plusEntries,
+                  plusLowerSlack, plusUpperSlack);
+    computeSlacks(col, step.minusHeaders, step.minusCoefs, step.minusEntries,
+                  minusLowerSlack, minusUpperSlack);
+
+    // non-basic propagation: if a generated row is non-basic (with nonzero
+    // dual), both its parents are forced non-basic. mark them so the greedy
+    // passes skip them. only force if the parent doesn't must-be-basic.
+    std::vector<HighsBool> forcedNonBasicPlus(numPlus, false);
+    std::vector<HighsBool> forcedNonBasicMinus(numMinus, false);
+    for (const auto& nr : step.newRows) {
+      // get indices of parent rows
+      HighsInt p = nr.plusParentIdx;
+      HighsInt m = nr.minusParentIdx;
+      // skip basic rows (zero dual) and degenerate non-basic rows (with zero
+      // dual)
+      if (p < 0 || m < 0 || std::abs(solution.row_dual[nr.row]) <= dual_tol)
+        continue;
+      // mark rows that do not have to be basic
+      if (!rowMustBeBasic(step.plusHeaders[p].row, plusLowerSlack[p],
+                          plusUpperSlack[p]))
+        forcedNonBasicPlus[p] = true;
+      if (!rowMustBeBasic(step.minusHeaders[m].row, minusLowerSlack[m],
+                          minusUpperSlack[m]))
+        forcedNonBasicMinus[m] = true;
+    }
+
+    // mark ranged rows (appearing in both plus and minus sets)
+    std::vector<HighsBool> isMinusRowRanged(numMinus, false);
+    HighsInt numRanged = 0;
+    for (HighsInt m = 0; m < numMinus; ++m)
+      for (HighsInt p = 0; p < numPlus; ++p)
+        if (step.minusHeaders[m].row == step.plusHeaders[p].row) {
+          isMinusRowRanged[m] = true;
+          numRanged++;
+          break;
+        }
+
+    // count number of basic new rows
+    HighsInt numNewRows = static_cast<HighsInt>(step.newRows.size());
+    HighsInt numBasicNewRows = 0;
+    for (const auto& nr : step.newRows)
+      if (basis.row_status[nr.row] == HighsBasisStatus::kBasic)
+        numBasicNewRows++;
+
+    // how many basic variables are needed?
+    HighsInt basicNeeded =
+        (numPlus + numMinus - numRanged - numNewRows) + numBasicNewRows;
+    HighsInt basicAssigned = 0;
+
+    // determine col status
+    bool colMustBeBasic =
+        solution.col_value[col] > step.header.colLower + tol &&
+        solution.col_value[col] < step.header.colUpper - tol;
+    bool colCanBeBasic =
+        colMustBeBasic || std::abs(solution.col_dual[col]) <= dual_tol;
+
+    // pass 1: assign all must-be-basic (col and rows)
+    if (colMustBeBasic) {
+      basis.col_status[col] = HighsBasisStatus::kBasic;
+      basicAssigned++;
+    }
+    for (HighsInt p = 0; p < numPlus; ++p) {
+      if (forcedNonBasicPlus[p]) continue;
+      if (rowMustBeBasic(step.plusHeaders[p].row, plusLowerSlack[p],
+                         plusUpperSlack[p]))
+        assignBasicRowStatus(step.plusHeaders[p].row, basicAssigned);
+    }
+    for (HighsInt m = 0; m < numMinus; ++m) {
+      if (isMinusRowRanged[m] || forcedNonBasicMinus[m]) continue;
+      if (rowMustBeBasic(step.minusHeaders[m].row, minusLowerSlack[m],
+                         minusUpperSlack[m]))
+        assignBasicRowStatus(step.minusHeaders[m].row, basicAssigned);
+    }
+
+    // pass 2: assign can-be-basic col (if not already assigned)
+    if (!colMustBeBasic) {
+      if (colCanBeBasic && basicAssigned < basicNeeded) {
+        basis.col_status[col] = HighsBasisStatus::kBasic;
+        basicAssigned++;
+      } else if (solution.col_value[col] <= step.header.colLower + tol) {
+        basis.col_status[col] = HighsBasisStatus::kLower;
+      } else {
+        basis.col_status[col] = HighsBasisStatus::kUpper;
+      }
+    }
+
+    // pass 3: assign can-be-basic rows, sorted by non-basic support count
+    // to reduce risk of rank deficiency in degenerate cases
+    std::vector<std::tuple<HighsInt, HighsInt, bool>> candidates;
+    for (HighsInt p = 0; p < numPlus; ++p)
+      collectCandidate(step.plusHeaders[p].row, forcedNonBasicPlus[p],
+                       plusLowerSlack[p], plusUpperSlack[p],
+                       step.plusEntries[p], p, false, candidates);
+    for (HighsInt m = 0; m < numMinus; ++m) {
+      if (isMinusRowRanged[m]) continue;
+      collectCandidate(step.minusHeaders[m].row, forcedNonBasicMinus[m],
+                       minusLowerSlack[m], minusUpperSlack[m],
+                       step.minusEntries[m], m, true, candidates);
+    }
+    // sort descending by non-basic support count
+    std::sort(candidates.begin(), candidates.end(),
+              [](const std::tuple<HighsInt, HighsInt, bool>& a,
+                 const std::tuple<HighsInt, HighsInt, bool>& b) {
+                return std::get<0>(a) > std::get<0>(b);
+              });
+    for (const auto& cand : candidates) {
+      HighsInt parentIndex = std::get<1>(cand);
+      if (std::get<2>(cand)) {
+        assignRowStatus(step.minusHeaders[parentIndex].row,
+                        minusLowerSlack[parentIndex],
+                        minusUpperSlack[parentIndex], basicAssigned,
+                        basicAssigned >= basicNeeded);
+      } else {
+        assignRowStatus(step.plusHeaders[parentIndex].row,
+                        plusLowerSlack[parentIndex],
+                        plusUpperSlack[parentIndex], basicAssigned,
+                        basicAssigned >= basicNeeded);
+      }
+    }
+  }
+}
+
 void HighsPostsolveStack::ZeroObjSingletonContinuousCol::undo(
     const HighsOptions& options, const std::vector<Nonzero>& rowValues,
     HighsSolution& solution, HighsBasis& basis) {
diff --git a/highs/presolve/HighsPostsolveStack.h b/highs/presolve/HighsPostsolveStack.h
index 1a9fdb0a7a..c236d3d03f 100644
--- a/highs/presolve/HighsPostsolveStack.h
+++ b/highs/presolve/HighsPostsolveStack.h
@@ -62,6 +62,55 @@ class HighsPostsolveStack {
     Nonzero() = default;
   };
 
+  template <typename RowStorageFormat>
+  struct FmeRowData {
+    HighsInt row;
+    double rowLower;
+    double rowUpper;
+    HighsMatrixSlice<RowStorageFormat> rowVec;
+  };
+
+  struct FmeRowHeader {
+    HighsInt row;
+    double rowLower;
+    double rowUpper;
+  };
+
+  struct FmeStepHeader {
+    double colLower;
+    double colUpper;
+    HighsInt col;
+    HighsInt numPlus;
+    HighsInt numMinus;
+  };
+
+  struct FmeDescendant {
+    HighsInt row;
+    double scaleFactor;
+  };
+
+  struct FmeNewRow {
+    HighsInt row;
+    HighsInt plusParentIdx;
+    HighsInt minusParentIdx;
+  };
+
+  struct FmeAncestryEntry {
+    HighsInt step;
+    HighsInt parentRowIndex;
+    double scale;
+    bool isMinus;
+  };
+
+  struct FmeBlockStep {
+    HighsInt col;
+    double colLower;
+    double colUpper;
+    HighsInt numPlus;
+    HighsInt numMinus;
+    std::vector<FmeNewRow> newRows;
+  };
+
   size_t debug_prev_numreductions = 0;
   double debug_prev_col_lower = 0;
   double debug_prev_col_upper = 0;
@@ -81,6 +130,17 @@ class HighsPostsolveStack {
     void transformToPresolvedSpace(std::vector<double>& primalSol) const;
   };
 
+  struct FourierMotzkinObjCol {
+    double offset;
+    HighsInt col;
+
+    void transformToPresolvedSpace(const std::vector<Nonzero>& costEntries,
+                                   std::vector<double>& primalSol) const;
+
+    void undo(const std::vector<Nonzero>& costEntries,
+              HighsSolution& solution) const;
+  };
+
   struct FreeColSubstitution {
     double rhs;
     double colCost;
@@ -273,8 +333,29 @@ class HighsPostsolveStack {
     kDuplicateColumn,
     kSlackColSubstitution,
     kZeroObjSingletonContinuousCol,
+    kFourierMotzkinBlock,
+    kFourierMotzkinObjCol,
   };
 
+  struct FmeStepData {
+    FmeStepHeader header;
+    std::vector<FmeRowHeader> plusHeaders;
+    std::vector<double> plusCoefs;
+    std::vector<std::vector<Nonzero>> plusEntries;
+    std::vector<std::vector<FmeDescendant>> plusDescendants;
+    std::vector<FmeRowHeader> minusHeaders;
+    std::vector<double> minusCoefs;
+    std::vector<std::vector<Nonzero>> minusEntries;
+    std::vector<std::vector<FmeDescendant>> minusDescendants;
+    std::vector<FmeNewRow> newRows;
+  };
+
+  static std::vector<FmeStepData> popFourierMotzkinBlock(HighsDataStack& stack);
+  static void undoFourierMotzkinBlock(const std::vector<FmeStepData>& steps,
+                                      const HighsOptions& options,
+                                      HighsSolution& solution,
+                                      HighsBasis& basis);
+
   HighsDataStack reductionValues;
   std::vector<std::pair<ReductionType, size_t>> reductions;
   std::vector<HighsInt> origColIndex;
@@ -345,6 +426,12 @@ class HighsPostsolveStack {
       case ReductionType::kZeroObjSingletonContinuousCol: {
         return "Zero obj singleton continuous col";
       }
+      case ReductionType::kFourierMotzkinBlock: {
+        return "Fourier-Motzkin block";
+      }
+      case ReductionType::kFourierMotzkinObjCol: {
+        return "Fourier-Motzkin obj col";
+      }
       default:
         return "Unknown";
     }
@@ -554,6 +641,121 @@ class HighsPostsolveStack {
     reductionAdded(ReductionType::kZeroObjSingletonContinuousCol);
   }
 
+  template <typename RowStorageFormat>
+  void fourierMotzkinBlockPushStep(
+      HighsInt col, const std::vector<FmeRowData<RowStorageFormat>>& plusRows,
+      const std::vector<FmeRowData<RowStorageFormat>>& minusRows) {
+    std::vector<FmeRowHeader> plusHeaders;
+    std::vector<double> plusCoefs;
+    plusHeaders.reserve(plusRows.size());
+    plusCoefs.reserve(plusRows.size());
+    for (const auto& rd : plusRows) {
+      std::vector<Nonzero> translated;
+      double coef = 0.0;
+      for (const HighsSliceNonzero& nz : rd.rowVec) {
+        if (nz.index() == col)
+          coef = nz.value();
+        else
+          translated.push_back({origColIndex[nz.index()], nz.value()});
+      }
+      reductionValues.push(translated);
+      plusCoefs.push_back(coef);
+      plusHeaders.push_back({origRowIndex[rd.row], rd.rowLower, rd.rowUpper});
+    }
+    reductionValues.push(plusCoefs);
+    reductionValues.push(plusHeaders);
+
+    std::vector<FmeRowHeader> minusHeaders;
+    std::vector<double> minusCoefs;
+    minusHeaders.reserve(minusRows.size());
+    minusCoefs.reserve(minusRows.size());
+    for (const auto& rd : minusRows) {
+      std::vector<Nonzero> translated;
+      double coef = 0.0;
+      for (const HighsSliceNonzero& nz : rd.rowVec) {
+        if (nz.index() == col)
+          coef = nz.value();
+        else
+          translated.push_back({origColIndex[nz.index()], nz.value()});
+      }
+      reductionValues.push(translated);
+      minusCoefs.push_back(coef);
+      minusHeaders.push_back({origRowIndex[rd.row], rd.rowLower, rd.rowUpper});
+    }
+    reductionValues.push(minusCoefs);
+    reductionValues.push(minusHeaders);
+  }
+
+  void fourierMotzkinBlockFinalise(
+      const std::vector<FmeBlockStep>& blockSteps,
+      const std::unordered_map<HighsInt, std::vector<FmeAncestryEntry>>&
+          rowAncestry) {
+    HighsInt numSteps = static_cast<HighsInt>(blockSteps.size());
+
+    std::vector<std::vector<std::vector<FmeDescendant>>> plusDescendantsAll(
+        numSteps);
+    std::vector<std::vector<std::vector<FmeDescendant>>> minusDescendantsAll(
+        numSteps);
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      plusDescendantsAll[s].resize(blockSteps[s].numPlus);
+      minusDescendantsAll[s].resize(blockSteps[s].numMinus);
+    }
+    for (const auto& entry : rowAncestry) {
+      HighsInt row = entry.first;
+      HighsInt origRow = origRowIndex[row];
+      for (const auto& a : entry.second) {
+        if (a.isMinus)
+          minusDescendantsAll[a.step][a.parentRowIndex].push_back(
+              {origRow, a.scale});
+        else
+          plusDescendantsAll[a.step][a.parentRowIndex].push_back(
+              {origRow, a.scale});
+      }
+    }
+
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      assert(static_cast<HighsInt>(plusDescendantsAll[s].size()) ==
+             blockSteps[s].numPlus);
+      for (HighsInt p = 0; p < blockSteps[s].numPlus; ++p)
+        reductionValues.push(plusDescendantsAll[s][p]);
+      assert(static_cast<HighsInt>(minusDescendantsAll[s].size()) ==
+             blockSteps[s].numMinus);
+      for (HighsInt m = 0; m < blockSteps[s].numMinus; ++m)
+        reductionValues.push(minusDescendantsAll[s][m]);
+    }
+
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      std::vector<FmeNewRow> translated;
+      translated.reserve(blockSteps[s].newRows.size());
+      for (const auto& nr : blockSteps[s].newRows)
+        translated.push_back(
+            {origRowIndex[nr.row], nr.plusParentIdx, nr.minusParentIdx});
+      reductionValues.push(translated);
+    }
+
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      FmeStepHeader header{blockSteps[s].colLower, blockSteps[s].colUpper,
+                           origColIndex[blockSteps[s].col],
+                           blockSteps[s].numPlus, blockSteps[s].numMinus};
+      reductionValues.push(header);
+    }
+
+    reductionValues.push(numSteps);
+    reductionAdded(ReductionType::kFourierMotzkinBlock);
+  }
+
+  void fourierMotzkinObjCol(HighsInt col, double offset,
+                            const std::vector<Nonzero>& costEntries) {
+    reductionValues.push(FourierMotzkinObjCol{offset, origColIndex[col]});
+    std::vector<Nonzero> translatedEntries;
+    translatedEntries.reserve(costEntries.size());
+    for (const Nonzero& entry : costEntries)
+      if (entry.index != col)
+        translatedEntries.emplace_back(origColIndex[entry.index], entry.value);
+    reductionValues.push(translatedEntries);
+    reductionAdded(ReductionType::kFourierMotzkinObjCol);
+  }
+
   template <typename ColStorageFormat>
   void doubletonEquation(HighsInt row, HighsInt colSubst, HighsInt col,
                          double coefSubst, double coef, double rhs,
@@ -1098,6 +1300,19 @@ class HighsPostsolveStack {
           reduction.undo(options, rowValues, solution, basis);
           break;
         }
+        case ReductionType::kFourierMotzkinBlock: {
+          auto steps = popFourierMotzkinBlock(reductionValues_);
+          undoFourierMotzkinBlock(steps, options, solution, basis);
+          break;
+        }
+        case ReductionType::kFourierMotzkinObjCol: {
+          std::vector<Nonzero> costEntries;
+          reductionValues_.pop(costEntries);
+          FourierMotzkinObjCol reduction;
+          reductionValues_.pop(reduction);
+          reduction.undo(costEntries, solution);
+          break;
+        }
         default:
           printf("Reduction case %d not handled\n",
                  int(reductions[i - 1].first));