diff --git a/highs/mip/HighsDebugSol.cpp b/highs/mip/HighsDebugSol.cpp
index f9d8432547..02dfa21242 100644
--- a/highs/mip/HighsDebugSol.cpp
+++ b/highs/mip/HighsDebugSol.cpp
@@ -68,10 +68,12 @@ void HighsDebugSol::activate() {
 
       HighsCDouble debugsolobj = 0.0;
       for (HighsInt i = 0; i != mipsolver->orig_model_->num_col_; ++i)
-        debugsolobj += mipsolver->orig_model_->col_cost_[i] *
-                       HighsCDouble(debugOrigSolution[i]);
+        debugsolobj +=
+            static_cast<HighsCDouble>(mipsolver->orig_model_->col_cost_[i]) *
+            debugOrigSolution[i];
 
-      debugSolObjective = double(debugsolobj + mipsolver->orig_model_->offset_);
+      debugSolObjective =
+          static_cast<double>(debugsolobj + mipsolver->orig_model_->offset_);
       debugSolActive = true;
       printf("debug sol active\n");
       registerDomain(mipsolver->mipdata_->getDomain());
@@ -101,9 +103,10 @@ void HighsDebugSol::shrink(const std::vector<HighsInt>& newColIndex) {
   debugSolution.resize(mipsolver->model_->num_col_);
   HighsCDouble debugsolobj = 0.0;
   for (HighsInt i = 0; i != mipsolver->model_->num_col_; ++i)
-    debugsolobj += mipsolver->model_->col_cost_[i] * debugSolution[i];
+    debugsolobj += static_cast<HighsCDouble>(mipsolver->model_->col_cost_[i]) *
+                   debugSolution[i];
 
-  debugSolObjective = double(debugsolobj);
+  debugSolObjective = static_cast<double>(debugsolobj);
 
   conflictingBounds.clear();
 }
diff --git a/highs/mip/HighsDomain.cpp b/highs/mip/HighsDomain.cpp
index f09d0f5b81..c8ee9539e1 100644
--- a/highs/mip/HighsDomain.cpp
+++ b/highs/mip/HighsDomain.cpp
@@ -78,12 +78,13 @@ HighsDomain::HighsDomain(HighsMipSolver& mipsolver) : mipsolver(&mipsolver) {
 }
 
 void HighsDomain::addCutpool(HighsCutPool& cutpool) {
-  HighsInt cutpoolindex = cutpoolpropagation.size();
+  HighsInt cutpoolindex = static_cast<HighsInt>(cutpoolpropagation.size());
   cutpoolpropagation.emplace_back(cutpoolindex, this, cutpool);
 }
 
 void HighsDomain::addConflictPool(HighsConflictPool& conflictPool) {
-  HighsInt conflictPoolIndex = conflictPoolPropagation.size();
+  HighsInt conflictPoolIndex =
+      static_cast<HighsInt>(conflictPoolPropagation.size());
   conflictPoolPropagation.emplace_back(conflictPoolIndex, this, conflictPool);
 }
 
@@ -182,7 +183,7 @@ void HighsDomain::ConflictPoolPropagation::conflictAdded(HighsInt conflict) {
   const std::vector<HighsDomainChange>& conflictEntries =
       conflictpool_->getConflictEntryVector();
 
-  if (HighsInt(conflictFlag_.size()) <= conflict) {
+  if (static_cast<HighsInt>(conflictFlag_.size()) <= conflict) {
     watchedLiterals_.resize(2 * conflict + 2);
     conflictFlag_.resize(conflict + 1);
   }
@@ -453,7 +454,7 @@ void HighsDomain::CutpoolPropagation::cutAdded(HighsInt cut, bool propagate) {
     const HighsInt* arindex = cutpool->getMatrix().getARindex();
     const double* arvalue = cutpool->getMatrix().getARvalue();
 
-    if (HighsInt(activitycuts_.size()) <= cut) {
+    if (static_cast<HighsInt>(activitycuts_.size()) <= cut) {
       activitycuts_.resize(cut + 1);
       activitycutsinf_.resize(cut + 1);
       propagatecutflags_.resize(cut + 1, 2);
@@ -469,7 +470,7 @@ void HighsDomain::CutpoolPropagation::cutAdded(HighsInt cut, bool propagate) {
     const HighsInt* arindex = cutpool->getMatrix().getARindex();
     const double* arvalue = cutpool->getMatrix().getARvalue();
 
-    if (HighsInt(activitycuts_.size()) <= cut) {
+    if (static_cast<HighsInt>(activitycuts_.size()) <= cut) {
       activitycuts_.resize(cut + 1);
       activitycutsinf_.resize(cut + 1);
       propagatecutflags_.resize(cut + 1, 2);
@@ -493,13 +494,14 @@ void HighsDomain::CutpoolPropagation::cutDeleted(
     return;
   }
 
-  if (cut < (HighsInt)propagatecutflags_.size()) propagatecutflags_[cut] |= 2;
+  if (cut < static_cast<HighsInt>(propagatecutflags_.size()))
+    propagatecutflags_[cut] |= 2;
 }
 
 void HighsDomain::CutpoolPropagation::markPropagateCut(HighsInt cut) {
   if (!propagatecutflags_[cut] &&
       (activitycutsinf_[cut] == 1 ||
-       (cutpool->getRhs()[cut] - double(activitycuts_[cut]) <=
+       (cutpool->getRhs()[cut] - activitycuts_[cut] <=
         capacityThreshold_[cut]))) {
     propagatecutinds_.push_back(cut);
     propagatecutflags_[cut] |= 1;
@@ -537,8 +539,7 @@ void HighsDomain::CutpoolPropagation::updateActivityLbChange(
           }
 
           if (activitycutsinf_[row] == 0 &&
-              activitycuts_[row] - cutpool->getRhs()[row] >
-                  domain->mipsolver->mipdata_->feastol) {
+              activitycuts_[row] - cutpool->getRhs()[row] > domain->feastol()) {
             domain->mipsolver->mipdata_->debugSolution.nodePruned(*domain);
             domain->infeasible_ = true;
             domain->infeasible_pos = domain->domchgstack_.size();
@@ -603,8 +604,7 @@ void HighsDomain::CutpoolPropagation::updateActivityUbChange(
           }
 
           if (activitycutsinf_[row] == 0 &&
-              activitycuts_[row] - cutpool->getRhs()[row] >
-                  domain->mipsolver->mipdata_->feastol) {
+              activitycuts_[row] - cutpool->getRhs()[row] > domain->feastol()) {
             domain->mipsolver->mipdata_->debugSolution.nodePruned(*domain);
             domain->infeasible_ = true;
             domain->infeasible_pos = domain->domchgstack_.size();
@@ -723,19 +723,21 @@ HighsDomain::ObjectivePropagation::ObjectivePropagation(HighsDomain* domain)
   }
 
   // add contribution of remaining objective nonzeros
-  const HighsInt numObjNz = objNonzeros.size();
+  const HighsInt numObjNz = static_cast<HighsInt>(objNonzeros.size());
   for (HighsInt i = partitionStarts[numPartitions]; i < numObjNz; ++i) {
     HighsInt col = objNonzeros[i];
     if (cost[col] > 0.0) {
       if (domain->col_lower_[col] == -kHighsInf)
         ++numInfObjLower;
       else
-        objectiveLower += domain->col_lower_[col] * cost[col];
+        objectiveLower +=
+            static_cast<HighsCDouble>(domain->col_lower_[col]) * cost[col];
     } else {
       if (domain->col_upper_[col] == kHighsInf)
         ++numInfObjLower;
       else
-        objectiveLower += domain->col_upper_[col] * cost[col];
+        objectiveLower +=
+            static_cast<HighsCDouble>(domain->col_upper_[col]) * cost[col];
     }
   }
 
@@ -748,7 +750,7 @@ void HighsDomain::ObjectivePropagation::getPropagationConstraint(
     HighsInt& len, double& rhs, HighsInt domchgCol) {
   const HighsInt numPartitions = objFunc->getNumCliquePartitions();
   inds = objFunc->getObjectiveNonzeros().data();
-  len = objFunc->getObjectiveNonzeros().size();
+  len = static_cast<HighsInt>(objFunc->getObjectiveNonzeros().size());
   if (numPartitions == 0) {
     vals = objFunc->getObjectiveValuesPacked().data();
     rhs = domain->mipsolver->mipdata_->upper_limit;
@@ -787,7 +789,7 @@ void HighsDomain::ObjectivePropagation::getPropagationConstraint(
   }
 
   vals = propagationConsBuffer.data();
-  rhs = double(tmpRhs);
+  rhs = static_cast<double>(tmpRhs);
 }
 
 void HighsDomain::ObjectivePropagation::recomputeCapacityThreshold() {
@@ -811,7 +813,7 @@ void HighsDomain::ObjectivePropagation::recomputeCapacityThreshold() {
   }
 
   const auto& objNonzeros = objFunc->getObjectiveNonzeros();
-  const HighsInt numObjNzs = objNonzeros.size();
+  const HighsInt numObjNzs = static_cast<HighsInt>(objNonzeros.size());
   for (HighsInt i = partitionStarts[numPartitions]; i < numObjNzs; ++i) {
     HighsInt col = objNonzeros[i];
 
@@ -845,12 +847,12 @@ void HighsDomain::ObjectivePropagation::updateActivityLbChange(
     if (oldbound == -kHighsInf)
       --numInfObjLower;
     else
-      objectiveLower -= oldbound * cost[col];
+      objectiveLower -= static_cast<HighsCDouble>(oldbound) * cost[col];
 
     if (newbound == -kHighsInf)
       ++numInfObjLower;
     else
-      objectiveLower += newbound * cost[col];
+      objectiveLower += static_cast<HighsCDouble>(newbound) * cost[col];
 
     debugCheckObjectiveLower();
 
@@ -965,12 +967,12 @@ void HighsDomain::ObjectivePropagation::updateActivityUbChange(
     if (oldbound == kHighsInf)
       --numInfObjLower;
     else
-      objectiveLower -= oldbound * cost[col];
+      objectiveLower -= static_cast<HighsCDouble>(oldbound) * cost[col];
 
     if (newbound == kHighsInf)
       ++numInfObjLower;
     else
-      objectiveLower += newbound * cost[col];
+      objectiveLower += static_cast<HighsCDouble>(newbound) * cost[col];
 
     debugCheckObjectiveLower();
 
@@ -1069,7 +1071,7 @@ bool HighsDomain::ObjectivePropagation::shouldBePropagated() const {
   if (domain->infeasible_) return false;
   double upperLimit = domain->mipsolver->mipdata_->upper_limit;
   if (upperLimit == kHighsInf) return false;
-  if (upperLimit - double(objectiveLower) > capacityThreshold) return false;
+  if (upperLimit - objectiveLower > capacityThreshold) return false;
 
   return true;
 }
@@ -1083,7 +1085,7 @@ void HighsDomain::ObjectivePropagation::debugCheckObjectiveLower() const {
   const auto& partitionStarts = objFunc->getCliquePartitionStarts();
   const auto& objNonzeros = objFunc->getObjectiveNonzeros();
 
-  const HighsInt numObjNzs = objNonzeros.size();
+  const HighsInt numObjNzs = static_cast<HighsInt>(objNonzeros.size());
   for (HighsInt i = 0; i < numPartitions; ++i) {
     HighsInt start = partitionStarts[i];
     HighsInt end = partitionStarts[i + 1];
@@ -1105,18 +1107,19 @@ void HighsDomain::ObjectivePropagation::debugCheckObjectiveLower() const {
     HighsInt col = objNonzeros[i];
     if (cost[col] > 0) {
       if (domain->col_lower_[col] > -kHighsInf)
-        lowerFromScratch += domain->col_lower_[col] * cost[col];
+        lowerFromScratch +=
+            static_cast<HighsCDouble>(domain->col_lower_[col]) * cost[col];
       else
         ++numInf;
     } else {
       if (domain->col_upper_[col] < kHighsInf)
-        lowerFromScratch += domain->col_upper_[col] * cost[col];
+        lowerFromScratch +=
+            static_cast<HighsCDouble>(domain->col_upper_[col]) * cost[col];
       else
         ++numInf;
     }
   }
-  assert(std::fabs(double(lowerFromScratch - objectiveLower)) <=
-         domain->feastol());
+  assert(abs(lowerFromScratch - objectiveLower) <= domain->feastol());
   assert(numInf == numInfObjLower);
 #endif
 }
@@ -1141,7 +1144,7 @@ void HighsDomain::ObjectivePropagation::propagate() {
     // Scan non-binary columns for infinite bound contribution until the one
     // column that contributes with an infinite bound is found which is the only
     // column that can be propagated
-    HighsInt numCol = objNonzeros.size();
+    HighsInt numCol = static_cast<HighsInt>(objNonzeros.size());
     for (HighsInt i = objFunc->getNumBinariesInObjective(); i < numCol; ++i) {
       HighsInt col = objNonzeros[i];
       if ((cost[col] > 0 && domain->col_lower_[col] != -kHighsInf) ||
@@ -1226,7 +1229,7 @@ void HighsDomain::ObjectivePropagation::propagate() {
 
       if (domain->infeasible_) break;
 
-      const HighsInt numObjNzs = objNonzeros.size();
+      const HighsInt numObjNzs = static_cast<HighsInt>(objNonzeros.size());
       for (HighsInt i = objFunc->getCliquePartitionStarts()[numPartitions];
            i < numObjNzs; ++i) {
         HighsInt col = objNonzeros[i];
@@ -1268,7 +1271,7 @@ void HighsDomain::computeMinActivity(HighsInt start, HighsInt end,
       HighsInt col = ARindex[j];
       double val = ARvalue[j];
 
-      assert(col < int(col_lower_.size()));
+      assert(col < static_cast<HighsInt>(col_lower_.size()));
 
       HighsInt tmp;
       double lb = getColLowerPos(col, infeasible_pos - 1, tmp);
@@ -1288,7 +1291,7 @@ void HighsDomain::computeMinActivity(HighsInt start, HighsInt end,
       HighsInt col = ARindex[j];
       double val = ARvalue[j];
 
-      assert(col < int(col_lower_.size()));
+      assert(col < static_cast<HighsInt>(col_lower_.size()));
 
       double contributionmin =
           activityContributionMin(val, col_lower_[col], col_upper_[col]);
@@ -1313,7 +1316,7 @@ void HighsDomain::computeMaxActivity(HighsInt start, HighsInt end,
       HighsInt col = ARindex[j];
       double val = ARvalue[j];
 
-      assert(col < int(col_lower_.size()));
+      assert(col < static_cast<HighsInt>(col_lower_.size()));
 
       HighsInt tmp;
       double lb = getColLowerPos(col, infeasible_pos - 1, tmp);
@@ -1332,7 +1335,7 @@ void HighsDomain::computeMaxActivity(HighsInt start, HighsInt end,
       HighsInt col = ARindex[j];
       double val = ARvalue[j];
 
-      assert(col < int(col_lower_.size()));
+      assert(col < static_cast<HighsInt>(col_lower_.size()));
 
       double contributionmin =
           activityContributionMax(val, col_lower_[col], col_upper_[col]);
@@ -1351,19 +1354,17 @@ double HighsDomain::adjustedUb(HighsInt col, HighsCDouble boundVal,
   double bound;
 
   if (mipsolver->isColIntegral(col)) {
-    bound = static_cast<double>(floor(boundVal + mipsolver->mipdata_->feastol));
+    bound = static_cast<double>(floor(boundVal + feastol()));
     accept = bound < col_upper_[col] &&
-             col_upper_[col] - bound >
-                 1000.0 * mipsolver->mipdata_->feastol * std::fabs(bound);
+             col_upper_[col] - bound > 1000.0 * feastol() * std::fabs(bound);
   } else {
-    if (std::fabs(double(boundVal) - col_lower_[col]) <=
-        mipsolver->mipdata_->epsilon)
+    if (abs(boundVal - col_lower_[col]) <= epsilon())
       bound = col_lower_[col];
     else
-      bound = double(boundVal);
+      bound = static_cast<double>(boundVal);
     if (col_upper_[col] == kHighsInf)
       accept = true;
-    else if (bound + 1000.0 * mipsolver->mipdata_->feastol < col_upper_[col]) {
+    else if (bound + 1000.0 * feastol() < col_upper_[col]) {
       double relativeImprove = col_upper_[col] - bound;
       if (col_lower_[col] != -kHighsInf)
         relativeImprove /= col_upper_[col] - col_lower_[col];
@@ -1383,19 +1384,17 @@ double HighsDomain::adjustedLb(HighsInt col, HighsCDouble boundVal,
   double bound;
 
   if (mipsolver->isColIntegral(col)) {
-    bound = static_cast<double>(ceil(boundVal - mipsolver->mipdata_->feastol));
+    bound = static_cast<double>(ceil(boundVal - feastol()));
     accept = bound > col_lower_[col] &&
-             bound - col_lower_[col] >
-                 1000.0 * mipsolver->mipdata_->feastol * std::fabs(bound);
+             bound - col_lower_[col] > 1000.0 * feastol() * std::fabs(bound);
   } else {
-    if (std::fabs(col_upper_[col] - double(boundVal)) <=
-        mipsolver->mipdata_->epsilon)
+    if (abs(col_upper_[col] - boundVal) <= epsilon())
       bound = col_upper_[col];
     else
-      bound = double(boundVal);
+      bound = static_cast<double>(boundVal);
     if (col_lower_[col] == -kHighsInf)
       accept = true;
-    else if (bound - 1000.0 * mipsolver->mipdata_->feastol > col_lower_[col]) {
+    else if (bound - 1000.0 * feastol() > col_lower_[col]) {
       double relativeImprove = bound - col_lower_[col];
       if (col_upper_[col] != kHighsInf)
         relativeImprove /= col_upper_[col] - col_lower_[col];
@@ -1416,7 +1415,7 @@ HighsInt HighsDomain::propagateRowUpper(const HighsInt* Rindex,
                                         const HighsCDouble& minactivity,
                                         HighsInt ninfmin,
                                         HighsDomainChange* boundchgs) const {
-  assert(std::isfinite(double(minactivity)));
+  assert(std::isfinite(static_cast<double>(minactivity)));
   if (ninfmin > 1) return 0;
   HighsInt numchgs = 0;
   for (HighsInt i = 0; i != Rlen; ++i) {
@@ -1432,8 +1431,7 @@ HighsInt HighsDomain::propagateRowUpper(const HighsInt* Rindex,
     }
 
     HighsCDouble boundVal = (Rupper - minresact) / Rvalue[i];
-    if (std::fabs(double(boundVal) * kHighsTiny) > mipsolver->mipdata_->feastol)
-      continue;
+    if (abs(boundVal) * kHighsTiny > feastol()) continue;
 
     if (Rvalue[i] > 0) {
       bool accept;
@@ -1460,7 +1458,7 @@ HighsInt HighsDomain::propagateRowLower(const HighsInt* Rindex,
                                         const HighsCDouble& maxactivity,
                                         HighsInt ninfmax,
                                         HighsDomainChange* boundchgs) const {
-  assert(std::isfinite(double(maxactivity)));
+  assert(std::isfinite(static_cast<double>(maxactivity)));
   if (ninfmax > 1) return 0;
   HighsInt numchgs = 0;
   for (HighsInt i = 0; i != Rlen; ++i) {
@@ -1476,8 +1474,7 @@ HighsInt HighsDomain::propagateRowLower(const HighsInt* Rindex,
     }
 
     HighsCDouble boundVal = (Rlower - maxresact) / Rvalue[i];
-    if (std::fabs(double(boundVal) * kHighsTiny) > mipsolver->mipdata_->feastol)
-      continue;
+    if (abs(boundVal) * kHighsTiny > feastol()) continue;
 
     if (Rvalue[i] < 0) {
       bool accept;
@@ -1501,14 +1498,12 @@ void HighsDomain::updateThresholdLbChange(HighsInt col, double newbound,
                                           double val, double& threshold) const {
   if (newbound != col_upper_[col]) {
     double thresholdNew =
-        std::fabs(val) * boundRange(col_upper_[col], newbound,
-                                    mipsolver->mipdata_->feastol,
-                                    variableType(col));
+        std::fabs(val) *
+        boundRange(col_upper_[col], newbound, feastol(), variableType(col));
 
     // the new threshold is now the maximum of the new threshold and the current
     // one
-    threshold =
-        std::max({threshold, thresholdNew, mipsolver->mipdata_->feastol});
+    threshold = std::max({threshold, thresholdNew, feastol()});
   }
 }
 
@@ -1516,14 +1511,12 @@ void HighsDomain::updateThresholdUbChange(HighsInt col, double newbound,
                                           double val, double& threshold) const {
   if (newbound != col_lower_[col]) {
     double thresholdNew =
-        std::fabs(val) * boundRange(newbound, col_lower_[col],
-                                    mipsolver->mipdata_->feastol,
-                                    variableType(col));
+        std::fabs(val) *
+        boundRange(newbound, col_lower_[col], feastol(), variableType(col));
 
     // the new threshold is now the maximum of the new threshold and the current
     // one
-    threshold =
-        std::max({threshold, thresholdNew, mipsolver->mipdata_->feastol});
+    threshold = std::max({threshold, thresholdNew, feastol()});
   }
 }
 
@@ -1556,8 +1549,8 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound,
             mipsolver->mipdata_->ARstart_[mip->a_matrix_.index_[i] + 1],
             mipsolver->mipdata_->ARindex_.data(),
             mipsolver->mipdata_->ARvalue_.data(), tmpinf, tmpminact);
-        assert(std::fabs(double(activitymin_[mip->a_matrix_.index_[i]] -
-                                tmpminact)) <= mipsolver->mipdata_->feastol);
+        assert(abs(activitymin_[mip->a_matrix_.index_[i]] - tmpminact) <=
+               feastol());
         assert(tmpinf == activitymininf_[mip->a_matrix_.index_[i]]);
       }
 #endif
@@ -1577,7 +1570,7 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound,
           activitymininf_[mip->a_matrix_.index_[i]] == 0 &&
           activitymin_[mip->a_matrix_.index_[i]] -
                   mip->row_upper_[mip->a_matrix_.index_[i]] >
-              mipsolver->mipdata_->feastol) {
+              feastol()) {
         mipsolver->mipdata_->debugSolution.nodePruned(*this);
         infeasible_ = true;
         infeasible_pos = domchgstack_.size();
@@ -1605,8 +1598,8 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound,
             mipsolver->mipdata_->ARstart_[mip->a_matrix_.index_[i] + 1],
             mipsolver->mipdata_->ARindex_.data(),
             mipsolver->mipdata_->ARvalue_.data(), tmpinf, tmpmaxact);
-        assert(std::fabs(double(activitymax_[mip->a_matrix_.index_[i]] -
-                                tmpmaxact)) <= mipsolver->mipdata_->feastol);
+        assert(abs(activitymax_[mip->a_matrix_.index_[i]] - tmpmaxact) <=
+               feastol());
         assert(tmpinf == activitymaxinf_[mip->a_matrix_.index_[i]]);
       }
 #endif
@@ -1626,7 +1619,7 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound,
           activitymaxinf_[mip->a_matrix_.index_[i]] == 0 &&
           mip->row_lower_[mip->a_matrix_.index_[i]] -
                   activitymax_[mip->a_matrix_.index_[i]] >
-              mipsolver->mipdata_->feastol) {
+              feastol()) {
         mipsolver->mipdata_->debugSolution.nodePruned(*this);
         infeasible_ = true;
         infeasible_pos = domchgstack_.size();
@@ -1723,8 +1716,8 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound,
             mipsolver->mipdata_->ARstart_[mip->a_matrix_.index_[i] + 1],
             mipsolver->mipdata_->ARindex_.data(),
             mipsolver->mipdata_->ARvalue_.data(), tmpinf, tmpmaxact);
-        assert(std::fabs(double(activitymax_[mip->a_matrix_.index_[i]] -
-                                tmpmaxact)) <= mipsolver->mipdata_->feastol);
+        assert(abs(activitymax_[mip->a_matrix_.index_[i]] - tmpmaxact) <=
+               feastol());
         assert(tmpinf == activitymaxinf_[mip->a_matrix_.index_[i]]);
       }
 #endif
@@ -1744,7 +1737,7 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound,
           activitymaxinf_[mip->a_matrix_.index_[i]] == 0 &&
           mip->row_lower_[mip->a_matrix_.index_[i]] -
                   activitymax_[mip->a_matrix_.index_[i]] >
-              mipsolver->mipdata_->feastol) {
+              feastol()) {
         mipsolver->mipdata_->debugSolution.nodePruned(*this);
         infeasible_ = true;
         infeasible_pos = domchgstack_.size();
@@ -1775,8 +1768,8 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound,
             mipsolver->mipdata_->ARstart_[mip->a_matrix_.index_[i] + 1],
             mipsolver->mipdata_->ARindex_.data(),
             mipsolver->mipdata_->ARvalue_.data(), tmpinf, tmpminact);
-        assert(std::fabs(double(activitymin_[mip->a_matrix_.index_[i]] -
-                                tmpminact)) <= mipsolver->mipdata_->feastol);
+        assert(abs(activitymin_[mip->a_matrix_.index_[i]] - tmpminact) <=
+               feastol());
         assert(tmpinf == activitymininf_[mip->a_matrix_.index_[i]]);
       }
 #endif
@@ -1796,7 +1789,7 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound,
           activitymininf_[mip->a_matrix_.index_[i]] == 0 &&
           activitymin_[mip->a_matrix_.index_[i]] -
                   mip->row_upper_[mip->a_matrix_.index_[i]] >
-              mipsolver->mipdata_->feastol) {
+              feastol()) {
         mipsolver->mipdata_->debugSolution.nodePruned(*this);
         infeasible_ = true;
         infeasible_pos = domchgstack_.size();
@@ -1923,9 +1916,10 @@ void HighsDomain::markPropagateCut(Reason reason) {
       break;
     default:
       assert(reason.type >= 0 &&
-             reason.type < HighsInt(cutpoolpropagation.size() +
-                                    conflictPoolPropagation.size()));
-      if (reason.type < (HighsInt)cutpoolpropagation.size())
+             reason.type <
+                 static_cast<HighsInt>(cutpoolpropagation.size() +
+                                       conflictPoolPropagation.size()));
+      if (reason.type < static_cast<HighsInt>(cutpoolpropagation.size()))
         cutpoolpropagation[reason.type].markPropagateCut(reason.index);
       else
         conflictPoolPropagation[reason.type - cutpoolpropagation.size()]
@@ -1935,20 +1929,20 @@ void HighsDomain::markPropagateCut(Reason reason) {
 
 void HighsDomain::markPropagate(HighsInt row) {
   if (!propagateflags_[row]) {
-    bool proplower = mipsolver->rowLower(row) != -kHighsInf &&
-                     (activitymininf_[row] != 0 ||
-                      activitymin_[row] < mipsolver->rowLower(row) -
-                                              mipsolver->mipdata_->feastol) &&
-                     (activitymaxinf_[row] == 1 ||
-                      (double(activitymax_[row]) - mipsolver->rowLower(row)) <=
-                          capacityThreshold_[row]);
-    bool propupper = mipsolver->rowUpper(row) != kHighsInf &&
-                     (activitymaxinf_[row] != 0 ||
-                      activitymax_[row] > mipsolver->rowUpper(row) +
-                                              mipsolver->mipdata_->feastol) &&
-                     (activitymininf_[row] == 1 ||
-                      (mipsolver->rowUpper(row) - double(activitymin_[row])) <=
-                          capacityThreshold_[row]);
+    bool proplower =
+        mipsolver->rowLower(row) != -kHighsInf &&
+        (activitymininf_[row] != 0 ||
+         activitymin_[row] < mipsolver->rowLower(row) - feastol()) &&
+        (activitymaxinf_[row] == 1 ||
+         activitymax_[row] - mipsolver->rowLower(row) <=
+             capacityThreshold_[row]);
+    bool propupper =
+        mipsolver->rowUpper(row) != kHighsInf &&
+        (activitymaxinf_[row] != 0 ||
+         activitymax_[row] > mipsolver->rowUpper(row) + feastol()) &&
+        (activitymininf_[row] == 1 ||
+         mipsolver->rowUpper(row) - activitymin_[row] <=
+             capacityThreshold_[row]);
 
     if (proplower || propupper) {
       propagateinds_.push_back(row);
@@ -2022,7 +2016,7 @@ double HighsDomain::doChangeBound(const HighsDomainChange& boundchg) {
 
 void HighsDomain::changeBound(HighsDomainChange boundchg, Reason reason) {
   assert(boundchg.column >= 0);
-  assert(boundchg.column < (HighsInt)col_upper_.size());
+  assert(boundchg.column < static_cast<HighsInt>(col_upper_.size()));
   // assert(infeasible_ == 0);
   // if (reason.type == Reason::kObjective) {
   //   if (!mipsolver->submip)
@@ -2038,8 +2032,7 @@ void HighsDomain::changeBound(HighsDomainChange boundchg, Reason reason) {
       boundchg.boundval = col_lower_[boundchg.column];
     }
     if (boundchg.boundval > col_upper_[boundchg.column]) {
-      if (boundchg.boundval - col_upper_[boundchg.column] >
-          mipsolver->mipdata_->feastol) {
+      if (boundchg.boundval - col_upper_[boundchg.column] > feastol()) {
         mipsolver->mipdata_->debugSolution.nodePruned(*this);
         if (!infeasible_) {
           infeasible_pos = domchgstack_.size();
@@ -2060,8 +2053,7 @@ void HighsDomain::changeBound(HighsDomainChange boundchg, Reason reason) {
       boundchg.boundval = col_upper_[boundchg.column];
     }
     if (boundchg.boundval < col_lower_[boundchg.column]) {
-      if (col_lower_[boundchg.column] - boundchg.boundval >
-          mipsolver->mipdata_->feastol) {
+      if (col_lower_[boundchg.column] - boundchg.boundval > feastol()) {
         mipsolver->mipdata_->debugSolution.nodePruned(*this);
         if (!infeasible_) {
           infeasible_pos = domchgstack_.size();
@@ -2084,7 +2076,7 @@ void HighsDomain::changeBound(HighsDomainChange boundchg, Reason reason) {
   if (reason.type == Reason::kBranching)
     branchPos_.push_back(domchgstack_.size());
 
-  assert(prevPos < (HighsInt)domchgstack_.size());
+  assert(prevPos < static_cast<HighsInt>(domchgstack_.size()));
 
   bool binary = isBinary(boundchg.column);
 
@@ -2106,8 +2098,7 @@ void HighsDomain::changeBound(HighsDomainChange boundchg, Reason reason) {
 
 bool HighsDomain::checkChangeBound(HighsBoundType boundtype, HighsInt col,
                                    HighsCDouble boundval, Reason reason) {
-  if (std::abs(static_cast<double>(boundval * kHighsTiny)) > feastol())
-    return false;
+  if (abs(boundval) * kHighsTiny > feastol()) return false;
   bool accept;
   double bound = boundtype == HighsBoundType::kLower
                      ? adjustedLb(col, boundval, accept)
@@ -2135,7 +2126,7 @@ void HighsDomain::setDomainChangeStack(
   domchgstack_.clear();
   domchgreason_.clear();
   branchPos_.clear();
-  HighsInt stacksize = domchgstack.size();
+  HighsInt stacksize = static_cast<HighsInt>(domchgstack.size());
   for (HighsInt k = 0; k != stacksize; ++k) {
     if (domchgstack[k].boundtype == HighsBoundType::kLower &&
         domchgstack[k].boundval <= col_lower_[domchgstack[k].column])
@@ -2169,7 +2160,7 @@ void HighsDomain::setDomainChangeStack(
   domchgstack_.clear();
   domchgreason_.clear();
   branchPos_.clear();
-  HighsInt stacksize = domchgstack.size();
+  HighsInt stacksize = static_cast<HighsInt>(domchgstack.size());
   HighsInt k = 0;
   for (HighsInt branchPos : branchingPositions) {
     for (; k < branchPos; ++k) {
@@ -2232,13 +2223,14 @@ void HighsDomain::setDomainChangeStack(
 }
 
 void HighsDomain::backtrackToGlobal() {
-  HighsInt k = HighsInt(domchgstack_.size()) - 1;
+  HighsInt k = static_cast<HighsInt>(domchgstack_.size()) - 1;
   bool old_infeasible = infeasible_;
   Reason old_reason = infeasible_reason;
 
-  if (infeasible_ && infeasible_pos == HighsInt(domchgstack_.size())) {
+  if (infeasible_ &&
+      infeasible_pos == static_cast<HighsInt>(domchgstack_.size())) {
     assert(old_infeasible);
-    assert(k == HighsInt(domchgstack_.size()) - 1);
+    assert(k == static_cast<HighsInt>(domchgstack_.size()) - 1);
     infeasible_ = false;
     infeasible_reason = Reason::unspecified();
   }
@@ -2267,7 +2259,7 @@ void HighsDomain::backtrackToGlobal() {
 
     if (infeasible_ && infeasible_pos == k) {
       assert(old_infeasible);
-      assert(k == HighsInt(domchgstack_.size()) - 1);
+      assert(k == static_cast<HighsInt>(domchgstack_.size()) - 1);
       infeasible_ = false;
       infeasible_reason = Reason::unspecified();
     }
@@ -2281,7 +2273,7 @@ void HighsDomain::backtrackToGlobal() {
     infeasible_ = false;
   }
 
-  HighsInt numreason = domchgreason_.size();
+  HighsInt numreason = static_cast<HighsInt>(domchgreason_.size());
   for (HighsInt i = k + 1; i < numreason; ++i)
     markPropagateCut(domchgreason_[i]);
 
@@ -2292,13 +2284,14 @@ void HighsDomain::backtrackToGlobal() {
 }
 
 HighsDomainChange HighsDomain::backtrack() {
-  HighsInt k = HighsInt(domchgstack_.size()) - 1;
+  HighsInt k = static_cast<HighsInt>(domchgstack_.size()) - 1;
   bool old_infeasible = infeasible_;
   Reason old_reason = infeasible_reason;
 
-  if (infeasible_ && infeasible_pos == HighsInt(domchgstack_.size())) {
+  if (infeasible_ &&
+      infeasible_pos == static_cast<HighsInt>(domchgstack_.size())) {
     assert(old_infeasible);
-    assert(k == HighsInt(domchgstack_.size()) - 1);
+    assert(k == static_cast<HighsInt>(domchgstack_.size()) - 1);
     infeasible_ = false;
     infeasible_reason = Reason::unspecified();
   }
@@ -2324,7 +2317,7 @@ HighsDomainChange HighsDomain::backtrack() {
 
     if (infeasible_ && infeasible_pos == k) {
       assert(old_infeasible);
-      assert(k == HighsInt(domchgstack_.size()) - 1);
+      assert(k == static_cast<HighsInt>(domchgstack_.size()) - 1);
       infeasible_ = false;
       infeasible_reason = Reason::unspecified();
     }
@@ -2343,7 +2336,7 @@ HighsDomainChange HighsDomain::backtrack() {
     infeasible_ = false;
   }
 
-  HighsInt numreason = domchgreason_.size();
+  HighsInt numreason = static_cast<HighsInt>(domchgreason_.size());
   for (HighsInt i = k + 1; i < numreason; ++i)
     markPropagateCut(domchgreason_[i]);
 
@@ -2396,7 +2389,8 @@ bool HighsDomain::propagate() {
   while (havePropagationRows()) {
     if (objProp_.isActive()) objProp_.propagate();
 
-    const HighsInt numConflictPools = conflictPoolPropagation.size();
+    const HighsInt numConflictPools =
+        static_cast<HighsInt>(conflictPoolPropagation.size());
     for (HighsInt conflictPool = 0; conflictPool < numConflictPools;
          ++conflictPool) {
       auto& conflictprop = conflictPoolPropagation[conflictPool];
@@ -2413,7 +2407,7 @@ bool HighsDomain::propagate() {
     if (!propagateinds_.empty()) {
       propagateinds.swap(propagateinds_);
 
-      HighsInt numproprows = propagateinds.size();
+      HighsInt numproprows = static_cast<HighsInt>(propagateinds.size());
       for (HighsInt i = 0; i != numproprows; ++i) {
         HighsInt row = propagateinds[i];
         propagateflags_[row] = 0;
@@ -2435,8 +2429,7 @@ bool HighsDomain::propagate() {
 
           if (mipsolver->rowUpper(i) != kHighsInf &&
               (activitymaxinf_[i] != 0 ||
-               activitymax_[i] >
-                   mipsolver->rowUpper(i) + mipsolver->mipdata_->feastol)) {
+               activitymax_[i] > mipsolver->rowUpper(i) + feastol())) {
             // computeMinActivity(start, end, mipsolver->ARstart_.data(),
             // mipsolver->ARvalue_.data(), activitymininf_[i],
             //           activitymin_[i]);
@@ -2450,8 +2443,7 @@ bool HighsDomain::propagate() {
 
           if (mipsolver->rowLower(i) != -kHighsInf &&
               (activitymininf_[i] != 0 ||
-               activitymin_[i] <
-                   mipsolver->rowLower(i) - mipsolver->mipdata_->feastol)) {
+               activitymin_[i] < mipsolver->rowLower(i) - feastol())) {
             // computeMaxActivity(start, end, mipsolver->ARstart_.data(),
             // mipsolver->ARvalue_.data(), activitymaxinf_[i],
             //           activitymax_[i]);
@@ -2497,13 +2489,13 @@ bool HighsDomain::propagate() {
       propagateinds.clear();
     }
 
-    const HighsInt numpools = cutpoolpropagation.size();
+    const HighsInt numpools = static_cast<HighsInt>(cutpoolpropagation.size());
     for (HighsInt cutpool = 0; cutpool != numpools; ++cutpool) {
       auto& cutpoolprop = cutpoolpropagation[cutpool];
       if (!cutpoolprop.propagatecutinds_.empty()) {
         propagateinds.swap(cutpoolprop.propagatecutinds_);
 
-        HighsInt numproprows = propagateinds.size();
+        HighsInt numproprows = static_cast<HighsInt>(propagateinds.size());
         for (HighsInt i = 0; i != numproprows; ++i) {
           HighsInt cut = propagateinds[i];
           cutpoolprop.propagatecutflags_[cut] &= 2;
@@ -2637,8 +2629,10 @@ void HighsDomain::conflictAnalyzeReconvergence(
 
   if (!conflictSet.explainBoundChangeLeq(
           conflictSet.reconvergenceFrontier,
-          ConflictSet::LocalDomChg{HighsInt(domchgstack_.size()), domchg},
-          proofinds, proofvals, prooflen, proofrhs, double(activitymin)))
+          ConflictSet::LocalDomChg{static_cast<HighsInt>(domchgstack_.size()),
+                                   domchg},
+          proofinds, proofvals, prooflen, proofrhs,
+          static_cast<double>(activitymin)))
     return;
 
   if (10 * conflictSet.resolvedDomainChanges.size() >
@@ -2649,7 +2643,7 @@ void HighsDomain::conflictAnalyzeReconvergence(
       conflictSet.resolvedDomainChanges.begin(),
       conflictSet.resolvedDomainChanges.end());
 
-  HighsInt depth = branchPos_.size();
+  HighsInt depth = static_cast<HighsInt>(branchPos_.size());
 
   while (depth > 0) {
     HighsInt branchPos = branchPos_[depth - 1];
@@ -2682,7 +2676,7 @@ void HighsDomain::tightenCoefficients(HighsInt* inds, double* vals,
   }
 
   HighsCDouble maxabscoef = maxactivity - rhs;
-  if (maxabscoef > mipsolver->mipdata_->feastol) {
+  if (maxabscoef > feastol()) {
     HighsCDouble upper = rhs;
     HighsInt tightened = 0;
     for (HighsInt i = 0; i != len; ++i) {
@@ -2690,12 +2684,12 @@ void HighsDomain::tightenCoefficients(HighsInt* inds, double* vals,
       if (vals[i] > maxabscoef) {
         HighsCDouble delta = vals[i] - maxabscoef;
         upper -= delta * col_upper_[inds[i]];
-        vals[i] = double(maxabscoef);
+        vals[i] = static_cast<double>(maxabscoef);
         ++tightened;
       } else if (vals[i] < -maxabscoef) {
         HighsCDouble delta = -vals[i] - maxabscoef;
         upper += delta * col_lower_[inds[i]];
-        vals[i] = -double(maxabscoef);
+        vals[i] = -static_cast<double>(maxabscoef);
         ++tightened;
       }
     }
@@ -2703,8 +2697,8 @@ void HighsDomain::tightenCoefficients(HighsInt* inds, double* vals,
     if (tightened != 0) {
       // printf("tightened %" HIGHSINT_FORMAT " coefficients, rhs changed from
       // %g to %g\n",
-      //       tightened, rhs, double(upper));
-      rhs = double(upper);
+      //       tightened, rhs, static_cast<double>(upper));
+      rhs = static_cast<double>(upper);
     }
   }
 }
@@ -2715,10 +2709,11 @@ double HighsDomain::getMinCutActivity(const HighsCutPool& cutpool,
     if (cutpoolprop.cutpool == &cutpool) {
       // assert((cutpoolprop.propagatecutflags_[cut] & 2) == 0);
 
-      return cut < (HighsInt)cutpoolprop.propagatecutflags_.size() &&
+      return cut < static_cast<HighsInt>(
+                       cutpoolprop.propagatecutflags_.size()) &&
                      (cutpoolprop.propagatecutflags_[cut] & 2) == 0 &&
                      cutpoolprop.activitycutsinf_[cut] == 0
-                 ? double(cutpoolprop.activitycuts_[cut])
+                 ? static_cast<double>(cutpoolprop.activitycuts_[cut])
                  : -kHighsInf;
     }
   }
@@ -2730,20 +2725,19 @@ bool HighsDomain::isFixing(const HighsDomainChange& domchg) const {
   double otherbound = domchg.boundtype == HighsBoundType::kUpper
                           ? col_lower_[domchg.column]
                           : col_upper_[domchg.column];
-  return std::fabs(domchg.boundval - otherbound) <=
-         mipsolver->mipdata_->epsilon;
+  return std::fabs(domchg.boundval - otherbound) <= epsilon();
 }
 
 HighsDomainChange HighsDomain::flip(const HighsDomainChange& domchg) const {
   if (domchg.boundtype == HighsBoundType::kLower) {
-    HighsDomainChange flipped{domchg.boundval - mipsolver->mipdata_->feastol,
-                              domchg.column, HighsBoundType::kUpper};
+    HighsDomainChange flipped{domchg.boundval - feastol(), domchg.column,
+                              HighsBoundType::kUpper};
     if (mipsolver->isColIntegral(domchg.column))
       flipped.boundval = std::floor(flipped.boundval);
     return flipped;
   } else {
-    HighsDomainChange flipped{domchg.boundval + mipsolver->mipdata_->feastol,
-                              domchg.column, HighsBoundType::kLower};
+    HighsDomainChange flipped{domchg.boundval + feastol(), domchg.column,
+                              HighsBoundType::kLower};
     if (mipsolver->isColIntegral(domchg.column))
       flipped.boundval = std::ceil(flipped.boundval);
     return flipped;
@@ -2752,6 +2746,8 @@ HighsDomainChange HighsDomain::flip(const HighsDomainChange& domchg) const {
 
 double HighsDomain::feastol() const { return mipsolver->mipdata_->feastol; }
 
+double HighsDomain::epsilon() const { return mipsolver->mipdata_->epsilon; }
+
 HighsDomain::ConflictSet::ConflictSet(HighsDomain& localdom_,
                                       const HighsDomain& globaldom_)
     : localdom(localdom_),
@@ -2800,7 +2796,8 @@ bool HighsDomain::ConflictSet::explainBoundChangeGeq(
       if (it != currentFrontier.end()) {
         cand.baseBound = it->domchg.boundval;
         if (cand.baseBound != globaldom.col_upper_[col])
-          M += vals[i] * (cand.baseBound - globaldom.col_upper_[col]);
+          M += vals[i] * (static_cast<HighsCDouble>(cand.baseBound) -
+                          globaldom.col_upper_[col]);
         if (cand.baseBound <= ub) continue;
       } else
         cand.baseBound = globaldom.col_upper_[col];
@@ -2817,7 +2814,8 @@ bool HighsDomain::ConflictSet::explainBoundChangeGeq(
       if (it != currentFrontier.end()) {
         cand.baseBound = it->domchg.boundval;
         if (cand.baseBound != globaldom.col_lower_[col])
-          M += vals[i] * (cand.baseBound - globaldom.col_lower_[col]);
+          M += vals[i] * (static_cast<HighsCDouble>(cand.baseBound) -
+                          globaldom.col_lower_[col]);
         if (cand.baseBound >= lb) continue;
       } else
         cand.baseBound = globaldom.col_lower_[col];
@@ -2847,16 +2845,16 @@ bool HighsDomain::ConflictSet::explainBoundChangeGeq(
     // in case of an integral variable the bound was rounded and can be
     // relaxed by 1-feastol. We use 1 - 10 * feastol for numerical safety.
     if (domchg.domchg.boundtype == HighsBoundType::kLower)
-      b0 -= (1.0 - 10 * localdom.mipsolver->mipdata_->feastol);
+      b0 -= (1.0 - 10 * localdom.feastol());
     else
-      b0 += (1.0 - 10 * localdom.mipsolver->mipdata_->feastol);
+      b0 += (1.0 - 10 * localdom.feastol());
   } else {
     // for a continuous variable we relax the bound by epsilon to
     // accommodate for tiny rounding errors
     if (domchg.domchg.boundtype == HighsBoundType::kLower)
-      b0 -= localdom.mipsolver->mipdata_->epsilon;
+      b0 -= localdom.epsilon();
     else
-      b0 += localdom.mipsolver->mipdata_->epsilon;
+      b0 += localdom.epsilon();
   }
 
   // now multiply the bound constraint with the coefficient value in the
@@ -2868,9 +2866,11 @@ bool HighsDomain::ConflictSet::explainBoundChangeGeq(
 
   // M is the current residual activity initially
   if (domchgVal < 0)
-    M -= domchgVal * globaldom.col_lower_[domchg.domchg.column];
+    M -= static_cast<HighsCDouble>(domchgVal) *
+         globaldom.col_lower_[domchg.domchg.column];
   else
-    M -= domchgVal * globaldom.col_upper_[domchg.domchg.column];
+    M -= static_cast<HighsCDouble>(domchgVal) *
+         globaldom.col_upper_[domchg.domchg.column];
 
   return resolveLinearGeq(M, Mupper, vals);
 }
@@ -2908,7 +2908,8 @@ bool HighsDomain::ConflictSet::explainBoundChangeLeq(
       if (it != currentFrontier.end()) {
         cand.baseBound = it->domchg.boundval;
         if (cand.baseBound != globaldom.col_lower_[col])
-          M += vals[i] * (cand.baseBound - globaldom.col_lower_[col]);
+          M += vals[i] * (static_cast<HighsCDouble>(cand.baseBound) -
+                          globaldom.col_lower_[col]);
         if (cand.baseBound >= lb) continue;
       } else
         cand.baseBound = globaldom.col_lower_[col];
@@ -2924,7 +2925,8 @@ bool HighsDomain::ConflictSet::explainBoundChangeLeq(
       if (it != currentFrontier.end()) {
         cand.baseBound = it->domchg.boundval;
         if (cand.baseBound != globaldom.col_upper_[col])
-          M += vals[i] * (cand.baseBound - globaldom.col_upper_[col]);
+          M += vals[i] * (static_cast<HighsCDouble>(cand.baseBound) -
+                          globaldom.col_upper_[col]);
         if (cand.baseBound <= ub) continue;
       } else
         cand.baseBound = globaldom.col_upper_[col];
@@ -2956,16 +2958,16 @@ bool HighsDomain::ConflictSet::explainBoundChangeLeq(
     // in case of an integral variable the bound was rounded and can be
     // relaxed by 1-feastol. We use 1 - 10 * feastol for numerical safety
     if (domchg.domchg.boundtype == HighsBoundType::kLower)
-      b0 -= (1.0 - 10 * localdom.mipsolver->mipdata_->feastol);
+      b0 -= (1.0 - 10 * localdom.feastol());
     else
-      b0 += (1.0 - 10 * localdom.mipsolver->mipdata_->feastol);
+      b0 += (1.0 - 10 * localdom.feastol());
   } else {
     // for a continuous variable we relax the bound by epsilon to
     // accommodate for tiny rounding errors
     if (domchg.domchg.boundtype == HighsBoundType::kLower)
-      b0 -= localdom.mipsolver->mipdata_->epsilon;
+      b0 -= localdom.epsilon();
     else
-      b0 += localdom.mipsolver->mipdata_->epsilon;
+      b0 += localdom.epsilon();
   }
 
   // now multiply the bound constraint with the coefficient value in the
@@ -2977,9 +2979,11 @@ bool HighsDomain::ConflictSet::explainBoundChangeLeq(
 
   // M is the global residual activity initially
   if (domchgVal < 0)
-    M -= domchgVal * globaldom.col_upper_[domchg.domchg.column];
+    M -= static_cast<HighsCDouble>(domchgVal) *
+         globaldom.col_upper_[domchg.domchg.column];
   else
-    M -= domchgVal * globaldom.col_lower_[domchg.domchg.column];
+    M -= static_cast<HighsCDouble>(domchgVal) *
+         globaldom.col_lower_[domchg.domchg.column];
 
   return resolveLinearLeq(M, Mlower, vals);
 }
@@ -2987,7 +2991,7 @@ bool HighsDomain::ConflictSet::explainBoundChangeLeq(
 bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
                                                 const double* vals) {
   resolvedDomainChanges.clear();
-  double covered = double(M - Mupper);
+  double covered = static_cast<double>(M - Mupper);
   if (covered > 0) {
     for (const auto& reasonDomchg : resolveBuffer) {
       LocalDomChg locdomchg;
@@ -2997,8 +3001,8 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
       resolvedDomainChanges.push_back(locdomchg);
       assert(resolvedDomainChanges.back().pos >= 0);
       assert(resolvedDomainChanges.back().pos <
-             (HighsInt)localdom.domchgstack_.size());
-      covered = double(M - Mupper);
+             static_cast<HighsInt>(localdom.domchgstack_.size()));
+      covered = static_cast<double>(M - Mupper);
       if (covered <= 0) break;
     }
 
@@ -3009,7 +3013,8 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
       // from the explanation
       // HighsInt numRelaxed = 0;
       // HighsInt numDropped = 0;
-      for (HighsInt k = resolvedDomainChanges.size() - 1; k >= 0; --k) {
+      for (HighsInt k = static_cast<HighsInt>(resolvedDomainChanges.size()) - 1;
+           k >= 0; --k) {
         ResolveCandidate& reasonDomchg = resolveBuffer[k];
         LocalDomChg& locdomchg = resolvedDomainChanges[k];
         HighsInt i = reasonDomchg.valuePos;
@@ -3017,8 +3022,8 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
         if (locdomchg.domchg.boundtype == HighsBoundType::kLower) {
           double lb = locdomchg.domchg.boundval;
           double glb = reasonDomchg.baseBound;
-          double relaxLb =
-              double(((Mupper - (M - reasonDomchg.delta)) / vals[i]) + glb);
+          double relaxLb = static_cast<double>(
+              ((Mupper - (M - reasonDomchg.delta)) / vals[i]) + glb);
           if (!localdom.mipsolver->isColContinuous(col))
             relaxLb = std::ceil(relaxLb);
 
@@ -3026,9 +3031,12 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
 
           locdomchg.domchg.boundval = relaxLb;
 
-          if (relaxLb - glb <= localdom.mipsolver->mipdata_->epsilon) {
+          // only drop if explanation stays valid without this entry
+          if (relaxLb - glb <= localdom.epsilon() &&
+              M - reasonDomchg.delta - Mupper <= localdom.feastol()) {
             // domain change can be fully removed from conflict
-            HighsInt last = resolvedDomainChanges.size() - 1;
+            HighsInt last =
+                static_cast<HighsInt>(resolvedDomainChanges.size()) - 1;
             std::swap(resolvedDomainChanges[last], resolvedDomainChanges[k]);
             resolvedDomainChanges.resize(last);
 
@@ -3039,17 +3047,17 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
               locdomchg.pos = localdom.prevboundval_[locdomchg.pos].second;
 
             // bound can be relaxed
-            M += vals[i] * (relaxLb - lb);
+            M += vals[i] * (static_cast<HighsCDouble>(relaxLb) - lb);
             // ++numRelaxed;
           }
 
-          covered = double(M - Mupper);
+          covered = static_cast<double>(M - Mupper);
           if (covered >= -localdom.feastol()) break;
         } else {
           double ub = locdomchg.domchg.boundval;
           double gub = reasonDomchg.baseBound;
-          double relaxUb =
-              double(((Mupper - (M - reasonDomchg.delta)) / vals[i]) + gub);
+          double relaxUb = static_cast<double>(
+              ((Mupper - (M - reasonDomchg.delta)) / vals[i]) + gub);
           if (!localdom.mipsolver->isColContinuous(col))
             relaxUb = std::floor(relaxUb);
 
@@ -3057,9 +3065,12 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
 
           locdomchg.domchg.boundval = relaxUb;
 
-          if (relaxUb - gub >= -localdom.mipsolver->mipdata_->epsilon) {
+          // only drop if explanation stays valid without this entry
+          if (relaxUb - gub >= -localdom.epsilon() &&
+              M - reasonDomchg.delta - Mupper <= localdom.feastol()) {
             // domain change can be fully removed from conflict
-            HighsInt last = resolvedDomainChanges.size() - 1;
+            HighsInt last =
+                static_cast<HighsInt>(resolvedDomainChanges.size()) - 1;
             std::swap(resolvedDomainChanges[last], resolvedDomainChanges[k]);
             resolvedDomainChanges.resize(last);
 
@@ -3070,11 +3081,11 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
             while (relaxUb >= localdom.prevboundval_[locdomchg.pos].first)
               locdomchg.pos = localdom.prevboundval_[locdomchg.pos].second;
 
-            M += vals[i] * (relaxUb - ub);
+            M += vals[i] * (static_cast<HighsCDouble>(relaxUb) - ub);
             // ++numRelaxed;
           }
 
-          covered = double(M - Mupper);
+          covered = static_cast<double>(M - Mupper);
           if (covered >= -localdom.feastol()) break;
         }
       }
@@ -3084,7 +3095,7 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
       //          (int)numRelaxed, (int)numDropped,
       //          (int)resolvedDomainChanges.size());
 
-      assert(covered <= localdom.mipsolver->mipdata_->feastol);
+      assert(covered <= localdom.feastol());
     }
   }
 
@@ -3094,7 +3105,7 @@ bool HighsDomain::ConflictSet::resolveLinearGeq(HighsCDouble M, double Mupper,
 bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
                                                 const double* vals) {
   resolvedDomainChanges.clear();
-  double covered = double(M - Mlower);
+  double covered = static_cast<double>(M - Mlower);
   if (covered < 0) {
     for (const auto& reasonDomchg : resolveBuffer) {
       LocalDomChg locdomchg;
@@ -3104,8 +3115,8 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
       resolvedDomainChanges.push_back(locdomchg);
       assert(resolvedDomainChanges.back().pos >= 0);
       assert(resolvedDomainChanges.back().pos <
-             (HighsInt)localdom.domchgstack_.size());
-      covered = double(M - Mlower);
+             static_cast<HighsInt>(localdom.domchgstack_.size()));
+      covered = static_cast<double>(M - Mlower);
       if (covered >= 0) break;
     }
 
@@ -3121,7 +3132,8 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
       // from the explanation
       // HighsInt numRelaxed = 0;
       // HighsInt numDropped = 0;
-      for (HighsInt k = resolvedDomainChanges.size() - 1; k >= 0; --k) {
+      for (HighsInt k = static_cast<HighsInt>(resolvedDomainChanges.size()) - 1;
+           k >= 0; --k) {
         ResolveCandidate& reasonDomchg = resolveBuffer[k];
         LocalDomChg& locdomchg = resolvedDomainChanges[k];
         HighsInt i = reasonDomchg.valuePos;
@@ -3129,8 +3141,8 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
         if (locdomchg.domchg.boundtype == HighsBoundType::kLower) {
           double lb = locdomchg.domchg.boundval;
           double glb = reasonDomchg.baseBound;
-          double relaxLb =
-              double(((Mlower - (M - reasonDomchg.delta)) / vals[i]) + glb);
+          double relaxLb = static_cast<double>(
+              ((Mlower - (M - reasonDomchg.delta)) / vals[i]) + glb);
           if (!localdom.mipsolver->isColContinuous(col))
             relaxLb = std::ceil(relaxLb);
 
@@ -3138,9 +3150,12 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
 
           locdomchg.domchg.boundval = relaxLb;
 
-          if (relaxLb - glb <= localdom.mipsolver->mipdata_->epsilon) {
+          // only drop if explanation stays valid without this entry
+          if (relaxLb - glb <= localdom.epsilon() &&
+              M - reasonDomchg.delta - Mlower >= -localdom.feastol()) {
             // domain change can be fully removed from conflict
-            HighsInt last = resolvedDomainChanges.size() - 1;
+            HighsInt last =
+                static_cast<HighsInt>(resolvedDomainChanges.size()) - 1;
             std::swap(resolvedDomainChanges[last], resolvedDomainChanges[k]);
             resolvedDomainChanges.resize(last);
 
@@ -3151,17 +3166,17 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
             while (relaxLb <= localdom.prevboundval_[locdomchg.pos].first)
               locdomchg.pos = localdom.prevboundval_[locdomchg.pos].second;
 
-            M += vals[i] * (relaxLb - lb);
+            M += vals[i] * (static_cast<HighsCDouble>(relaxLb) - lb);
             // ++numRelaxed;
           }
 
-          covered = double(M - Mlower);
+          covered = static_cast<double>(M - Mlower);
           if (covered <= localdom.feastol()) break;
         } else {
           double ub = locdomchg.domchg.boundval;
           double gub = reasonDomchg.baseBound;
-          double relaxUb =
-              double(((Mlower - (M - reasonDomchg.delta)) / vals[i]) + gub);
+          double relaxUb = static_cast<double>(
+              ((Mlower - (M - reasonDomchg.delta)) / vals[i]) + gub);
           if (!localdom.mipsolver->isColContinuous(col))
             relaxUb = std::floor(relaxUb);
 
@@ -3169,9 +3184,12 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
 
           locdomchg.domchg.boundval = relaxUb;
 
-          if (relaxUb - gub >= -localdom.mipsolver->mipdata_->epsilon) {
+          // only drop if explanation stays valid without this entry
+          if (relaxUb - gub >= -localdom.epsilon() &&
+              M - reasonDomchg.delta - Mlower >= -localdom.feastol()) {
             // domain change can be fully removed from conflict
-            HighsInt last = resolvedDomainChanges.size() - 1;
+            HighsInt last =
+                static_cast<HighsInt>(resolvedDomainChanges.size()) - 1;
             std::swap(resolvedDomainChanges[last], resolvedDomainChanges[k]);
             resolvedDomainChanges.resize(last);
 
@@ -3182,11 +3200,11 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
             while (relaxUb >= localdom.prevboundval_[locdomchg.pos].first)
               locdomchg.pos = localdom.prevboundval_[locdomchg.pos].second;
 
-            M += vals[i] * (relaxUb - ub);
+            M += vals[i] * (static_cast<HighsCDouble>(relaxUb) - ub);
             // ++numRelaxed;
           }
 
-          covered = double(M - Mlower);
+          covered = static_cast<double>(M - Mlower);
           if (covered <= localdom.feastol()) break;
         }
       }
@@ -3196,7 +3214,7 @@ bool HighsDomain::ConflictSet::resolveLinearLeq(HighsCDouble M, double Mlower,
       //          (int)numRelaxed, (int)numDropped,
       //          (int)resolvedDomainChanges.size());
 
-      assert(covered >= -localdom.mipsolver->mipdata_->feastol);
+      assert(covered >= -localdom.feastol());
     }
   }
   return true;
@@ -3221,12 +3239,12 @@ bool HighsDomain::ConflictSet::explainInfeasibility() {
         double ub =
             localdom.getColUpperPos(col, conflictingBoundPos, otherBoundPos);
         assert(localdom.domchgstack_[conflictingBoundPos].boundval - ub >
-               +localdom.mipsolver->mipdata_->feastol);
+               +localdom.feastol());
       } else {
         double lb =
             localdom.getColLowerPos(col, conflictingBoundPos, otherBoundPos);
         assert(localdom.domchgstack_[conflictingBoundPos].boundval - lb <
-               -localdom.mipsolver->mipdata_->feastol);
+               -localdom.feastol());
       }
       if (otherBoundPos != -1)
         resolvedDomainChanges.push_back(
@@ -3281,16 +3299,17 @@ bool HighsDomain::ConflictSet::explainInfeasibility() {
       // heuristic.
       if (ninfmin > 0) return false;
 
-      return explainInfeasibilityLeq(inds, vals, len, rhs, double(minAct));
+      return explainInfeasibilityLeq(inds, vals, len, rhs,
+                                     static_cast<double>(minAct));
     }
     default:
       assert(localdom.infeasible_reason.type >= 0);
       assert(localdom.infeasible_reason.type <
-             HighsInt(localdom.cutpoolpropagation.size() +
-                      localdom.conflictPoolPropagation.size()));
+             static_cast<HighsInt>(localdom.cutpoolpropagation.size() +
+                                   localdom.conflictPoolPropagation.size()));
 
       if (localdom.infeasible_reason.type <
-          (HighsInt)localdom.cutpoolpropagation.size()) {
+          static_cast<HighsInt>(localdom.cutpoolpropagation.size())) {
         HighsInt cutpoolIndex = localdom.infeasible_reason.type;
         HighsInt cutIndex = localdom.infeasible_reason.index;
 
@@ -3413,8 +3432,7 @@ bool HighsDomain::ConflictSet::explainInfeasibilityGeq(const HighsInt* inds,
   pdqsort(resolveBuffer.begin(), resolveBuffer.end());
 
   // compute the lower bound of M that is necessary
-  double Mupper = rhs - std::max(10.0, std::fabs(rhs)) *
-                            localdom.mipsolver->mipdata_->feastol;
+  double Mupper = rhs - std::max(10.0, std::fabs(rhs)) * localdom.feastol();
 
   assert(reasonSideFrontier.empty());
   return resolveLinearGeq(maxAct, Mupper, vals);
@@ -3460,8 +3478,7 @@ bool HighsDomain::ConflictSet::explainInfeasibilityLeq(const HighsInt* inds,
   pdqsort(resolveBuffer.begin(), resolveBuffer.end());
 
   // compute the lower bound of M that is necessary
-  double Mlower = rhs + std::max(10.0, std::fabs(rhs)) *
-                            localdom.mipsolver->mipdata_->feastol;
+  double Mlower = rhs + std::max(10.0, std::fabs(rhs)) * localdom.feastol();
 
   return resolveLinearLeq(minAct, Mlower, vals);
 }
@@ -3481,13 +3498,13 @@ bool HighsDomain::ConflictSet::explainBoundChange(
       if (val) {
         assert(localdom.colLowerPos_[col] >= 0);
         assert(localdom.colLowerPos_[col] <
-               (HighsInt)localdom.domchgstack_.size());
+               static_cast<HighsInt>(localdom.domchgstack_.size()));
 
         localdom.getColLowerPos(col, domchg.pos, boundPos);
       } else {
         assert(localdom.colUpperPos_[col] >= 0);
         assert(localdom.colUpperPos_[col] <
-               (HighsInt)localdom.domchgstack_.size());
+               static_cast<HighsInt>(localdom.domchgstack_.size()));
 
         localdom.getColUpperPos(col, domchg.pos, boundPos);
       }
@@ -3546,15 +3563,15 @@ bool HighsDomain::ConflictSet::explainBoundChange(
       if (ninfmin == 1) return false;
 
       return explainBoundChangeLeq(currentFrontier, domchg, inds, vals, len,
-                                   rhs, double(minAct));
+                                   rhs, static_cast<double>(minAct));
     }
     default:
       assert(localdom.domchgreason_[domchg.pos].type >= 0);
       assert(localdom.domchgreason_[domchg.pos].type <
-             (HighsInt)(localdom.cutpoolpropagation.size() +
-                        localdom.conflictPoolPropagation.size()));
+             static_cast<HighsInt>(localdom.cutpoolpropagation.size() +
+                                   localdom.conflictPoolPropagation.size()));
       if (localdom.domchgreason_[domchg.pos].type <
-          (HighsInt)localdom.cutpoolpropagation.size()) {
+          static_cast<HighsInt>(localdom.cutpoolpropagation.size())) {
         HighsInt cutpoolIndex = localdom.domchgreason_[domchg.pos].type;
         HighsInt cutIndex = localdom.domchgreason_[domchg.pos].index;
 
@@ -3687,7 +3704,7 @@ HighsInt HighsDomain::ConflictSet::queueSize() const {
 
 bool HighsDomain::ConflictSet::resolvable(HighsInt domChgPos) const {
   assert(domChgPos >= 0);
-  assert(domChgPos < (HighsInt)localdom.domchgreason_.size());
+  assert(domChgPos < static_cast<HighsInt>(localdom.domchgreason_.size()));
   // printf("domchgPos: %d\n", domChgPos);
   // printf("stacksize: %ld\n", localdom.domchgreason_.size());
   switch (localdom.domchgreason_[domChgPos].type) {
@@ -3709,7 +3726,7 @@ HighsInt HighsDomain::ConflictSet::resolveDepth(std::set<LocalDomChg>& frontier,
   LocalDomChg startPos =
       LocalDomChg{depthLevel == 0 ? 0 : localdom.branchPos_[depthLevel - 1] + 1,
                   HighsDomainChange()};
-  while (depthLevel < (HighsInt)localdom.branchPos_.size()) {
+  while (depthLevel < static_cast<HighsInt>(localdom.branchPos_.size())) {
     HighsInt branchPos = localdom.branchPos_[depthLevel];
     if (localdom.domchgstack_[branchPos].boundval !=
         localdom.prevboundval_[branchPos].first)
@@ -3719,7 +3736,7 @@ HighsInt HighsDomain::ConflictSet::resolveDepth(std::set<LocalDomChg>& frontier,
   }
 
   auto iterEnd =
-      depthLevel == (HighsInt)localdom.branchPos_.size()
+      depthLevel == static_cast<HighsInt>(localdom.branchPos_.size())
           ? frontier.end()
           : frontier.upper_bound(LocalDomChg{localdom.branchPos_[depthLevel],
                                              HighsDomainChange()});
@@ -3787,7 +3804,8 @@ HighsInt HighsDomain::ConflictSet::computeCuts(HighsInt depthLevel,
                                                HighsPseudocost& pseudocost) {
   HighsInt numResolved = resolveDepth(
       reasonSideFrontier, depthLevel, 1, pseudocost,
-      depthLevel == (HighsInt)localdom.branchPos_.size() ? 1 : 0, true);
+      depthLevel == static_cast<HighsInt>(localdom.branchPos_.size()) ? 1 : 0,
+      true);
   if (numResolved == -1) return -1;
   HighsInt numConflicts = 0;
   if (numResolved > 0) {
@@ -3851,7 +3869,7 @@ void HighsDomain::ConflictSet::conflictAnalysis(HighsConflictPool& conflictPool,
       reasonSideFrontier, localdom.domchgstack_);
 
   HighsInt numConflicts = 0;
-  HighsInt lastDepth = localdom.branchPos_.size();
+  HighsInt lastDepth = static_cast<HighsInt>(localdom.branchPos_.size());
   // printf("start conflict analysis\n");
   HighsInt currDepth;
   for (currDepth = lastDepth; currDepth >= 0; --currDepth) {
@@ -3904,7 +3922,7 @@ void HighsDomain::ConflictSet::conflictAnalysis(const HighsInt* proofinds,
   if (ninfmin != 0) return;
 
   if (!explainInfeasibilityLeq(proofinds, proofvals, prooflen, proofrhs,
-                               double(activitymin)))
+                               static_cast<double>(activitymin)))
     return;
 
   pseudocost.increaseConflictWeight();
@@ -3928,7 +3946,7 @@ void HighsDomain::ConflictSet::conflictAnalysis(const HighsInt* proofinds,
       reasonSideFrontier, localdom.domchgstack_);
 
   HighsInt numConflicts = 0;
-  HighsInt lastDepth = localdom.branchPos_.size();
+  HighsInt lastDepth = static_cast<HighsInt>(localdom.branchPos_.size());
   HighsInt currDepth;
   for (currDepth = lastDepth; currDepth >= 0; --currDepth) {
     if (currDepth > 0) {
diff --git a/highs/mip/HighsDomain.h b/highs/mip/HighsDomain.h
index e69ab43c8c..95c742eb6e 100644
--- a/highs/mip/HighsDomain.h
+++ b/highs/mip/HighsDomain.h
@@ -537,7 +537,7 @@ class HighsDomain {
 
   double getObjectiveLowerBound() const {
     if (objProp_.isActive() && objProp_.numInfObjLower == 0)
-      return double(objProp_.objectiveLower);
+      return static_cast<double>(objProp_.objectiveLower);
 
     return -kHighsInf;
   }
@@ -561,7 +561,7 @@ class HighsDomain {
     std::vector<HighsDomainChange> reducedstack;
     reducedstack.reserve(domchgstack_.size());
     branchingPositions.reserve(branchPos_.size());
-    for (HighsInt i = 0; i < (HighsInt)domchgstack_.size(); ++i) {
+    for (HighsInt i = 0; i < static_cast<HighsInt>(domchgstack_.size()); ++i) {
       // keep only the tightest bound change for each variable
       if ((domchgstack_[i].boundtype == HighsBoundType::kLower &&
            colLowerPos_[domchgstack_[i].column] != i) ||
@@ -620,11 +620,13 @@ class HighsDomain {
                            double& rhs) const;
 
   double getMinActivity(HighsInt row) const {
-    return activitymininf_[row] == 0 ? double(activitymin_[row]) : -kHighsInf;
+    return activitymininf_[row] == 0 ? static_cast<double>(activitymin_[row])
+                                     : -kHighsInf;
   }
 
   double getMaxActivity(HighsInt row) const {
-    return activitymaxinf_[row] == 0 ? double(activitymax_[row]) : kHighsInf;
+    return activitymaxinf_[row] == 0 ? static_cast<double>(activitymax_[row])
+                                     : kHighsInf;
   }
 
   double getMinCutActivity(const HighsCutPool& cutpool, HighsInt cut) const;
@@ -654,6 +656,8 @@ class HighsDomain {
 
   double feastol() const;
 
+  double epsilon() const;
+
   HighsInt numModelNonzeros() const { return mipsolver->numNonzero(); }
 
   bool inSubmip() const { return mipsolver->submip; }
diff --git a/highs/mip/HighsLpRelaxation.cpp b/highs/mip/HighsLpRelaxation.cpp
index 77b2b0f98f..843b27849b 100644
--- a/highs/mip/HighsLpRelaxation.cpp
+++ b/highs/mip/HighsLpRelaxation.cpp
@@ -446,7 +446,7 @@ double HighsLpRelaxation::computeBestEstimate(const HighsPseudocost& ps) const {
     estimate += double(increase);
   }
 
-  return double(estimate);
+  return static_cast<double>(estimate);
 }
 
 double HighsLpRelaxation::computeLPDegneracy(
@@ -1612,7 +1612,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
           }
 
           for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
-            objsum += roundsol[i] * mipsolver.colCost(i);
+            objsum +=
+                static_cast<HighsCDouble>(roundsol[i]) * mipsolver.colCost(i);
 
           if (!mipsolver.mipdata_->parallelLockActive() || !worker_) {
             mipsolver.mipdata_->addIncumbent(
@@ -1625,7 +1626,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
         }
 
         for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
-          objsum += sol.col_value[i] * mipsolver.colCost(i);
+          objsum += static_cast<HighsCDouble>(sol.col_value[i]) *
+                    mipsolver.colCost(i);
 
         if (fractionalints.empty() && !unscaledPrimalFeasible(status)) {
           std::vector<double> fixSol = sol.col_value;
@@ -1647,7 +1649,7 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
           }
         }
 
-        objective = double(objsum);
+        objective = static_cast<double>(objsum);
         break;
       }
       case Status::kInfeasible:
diff --git a/highs/mip/HighsMipSolver.cpp b/highs/mip/HighsMipSolver.cpp
index 1742cc730d..2044074f48 100644
--- a/highs/mip/HighsMipSolver.cpp
+++ b/highs/mip/HighsMipSolver.cpp
@@ -1365,7 +1365,7 @@ bool HighsMipSolver::solutionFeasible(const HighsLp* lp,
     assert(col_value.size() == static_cast<size_t>(lp->num_col_));
   for (HighsInt i = 0; i != lp->num_col_; ++i) {
     const double value = col_value[i];
-    obj += lp->col_cost_[i] * value;
+    obj += static_cast<HighsCDouble>(lp->col_cost_[i]) * value;
 
     if (lp->integrality_[i] == HighsVarType::kInteger) {
       integrality_violation =
diff --git a/highs/mip/HighsMipSolverData.cpp b/highs/mip/HighsMipSolverData.cpp
index cabbdb8070..2b7bd3184e 100644
--- a/highs/mip/HighsMipSolverData.cpp
+++ b/highs/mip/HighsMipSolverData.cpp
@@ -226,7 +226,7 @@ bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
     if (mipsolver.isColInteger(i) && fractionality(solution[i]) > feastol)
       return false;
 
-    obj += mipsolver.colCost(i) * solution[i];
+    obj += static_cast<HighsCDouble>(mipsolver.colCost(i)) * solution[i];
   }
 
   for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
@@ -389,8 +389,9 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() {
 
     HighsCDouble cdouble_obj = 0.0;
     for (HighsInt iCol = 0; iCol < mipsolver.numCol(); iCol++)
-      cdouble_obj += mipsolver.colCost(iCol) * solution[iCol];
-    double obj = double(cdouble_obj);
+      cdouble_obj +=
+          static_cast<HighsCDouble>(mipsolver.colCost(iCol)) * solution[iCol];
+    double obj = static_cast<double>(cdouble_obj);
     const double save_upper_bound = upper_bound;
     const bool new_incumbent =
         addIncumbent(solution, obj, heuristic_source[try_heuristic]);
@@ -1145,8 +1146,8 @@ void HighsMipSolverData::runSetup() {
     debugSolution.debugSolObjective = 0;
     HighsCDouble debugsolobj = 0.0;
     for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
-      debugsolobj +=
-          mipsolver.colCost(i) * HighsCDouble(debugSolution.debugSolution[i]);
+      debugsolobj += static_cast<HighsCDouble>(mipsolver.colCost(i)) *
+                     debugSolution.debugSolution[i];
     debugSolution.debugSolObjective = static_cast<double>(debugsolobj);
     debugSolution.registerDomain(getDomain());
     assert(checkSolution(debugSolution.debugSolution));
diff --git a/highs/mip/HighsMipWorker.cpp b/highs/mip/HighsMipWorker.cpp
index 1fcba77573..7c7f2cd6a1 100644
--- a/highs/mip/HighsMipWorker.cpp
+++ b/highs/mip/HighsMipWorker.cpp
@@ -127,7 +127,7 @@ bool HighsMipWorker::trySolution(const std::vector<double>& solution,
         fractionality(solution[i]) > mipdata_.feastol)
       return false;
 
-    obj += mipsolver_.colCost(i) * solution[i];
+    obj += static_cast<HighsCDouble>(mipsolver_.colCost(i)) * solution[i];
   }
 
   for (HighsInt i = 0; i != mipsolver_.model_->num_row_; ++i) {
diff --git a/highs/mip/HighsSearch.cpp b/highs/mip/HighsSearch.cpp
index 5149c2f6f2..2ab0050823 100644
--- a/highs/mip/HighsSearch.cpp
+++ b/highs/mip/HighsSearch.cpp
@@ -44,7 +44,7 @@ double HighsSearch::checkSol(const std::vector<double>& sol,
   HighsCDouble objval = 0.0;
   integerfeasible = true;
   for (HighsInt i = 0; i != mipsolver.numCol(); ++i) {
-    objval += sol[i] * mipsolver.colCost(i);
+    objval += static_cast<HighsCDouble>(sol[i]) * mipsolver.colCost(i);
     assert(std::isfinite(sol[i]));
 
     if (!integerfeasible || !mipsolver.isColInteger(i)) continue;
diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp
index ac5b83a74f..954e432d2f 100644
--- a/highs/presolve/HPresolve.cpp
+++ b/highs/presolve/HPresolve.cpp
@@ -2952,7 +2952,8 @@ void HPresolve::substitute(HighsInt row, HighsInt col, double rhs) {
 
   // substitute column in the objective function
   if (model->col_cost_[col] != 0.0) {
-    HighsCDouble objscale = model->col_cost_[col] * substrowscale;
+    HighsCDouble objscale =
+        static_cast<HighsCDouble>(model->col_cost_[col]) * substrowscale;
     model->offset_ = static_cast<double>(model->offset_ - objscale * rhs);
     assert(std::isfinite(model->offset_));
     for (HighsInt rowiter : rowpositions) {
@@ -2967,7 +2968,7 @@ void HPresolve::substitute(HighsInt row, HighsInt col, double rhs) {
     }
     assert(std::abs(model->col_cost_[col]) <=
            std::max(options->dual_feasibility_tolerance,
-                    kHighsTiny * std::abs(static_cast<double>(objscale))));
+                    static_cast<double>(kHighsTiny * abs(objscale))));
     model->col_cost_[col] = 0.0;
   }
 
@@ -7104,10 +7105,12 @@ HPresolve::Result HPresolve::strengthenInequalities(
       double ub = model->col_upper_[col] - model->col_lower_[col];
       if (weight > 0) {
         comp = 1;
-        maxviolation += model->col_upper_[col] * weight;
+        maxviolation +=
+            static_cast<HighsCDouble>(model->col_upper_[col]) * weight;
       } else {
         comp = -1;
-        maxviolation += model->col_lower_[col] * weight;
+        maxviolation +=
+            static_cast<HighsCDouble>(model->col_lower_[col]) * weight;
         weight = -weight;
       }
 
@@ -7227,16 +7230,15 @@ HPresolve::Result HPresolve::strengthenInequalities(
                               HighsInt direction) {
       for (HighsInt i : indices) {
         assert(Arow[positions[i]] == row);
-        double coefdelta =
-            direction * static_cast<double>(reducedcost[i] - maxviolation);
+        HighsCDouble coefdelta = direction * (reducedcost[i] - maxviolation);
         HighsInt col = Acol[positions[i]];
 
         if (complementation[i] == -1) {
           rhs += coefdelta * model->col_lower_[col];
-          addToMatrix(row, col, coefdelta);
+          addToMatrix(row, col, static_cast<double>(coefdelta));
         } else {
           rhs -= coefdelta * model->col_upper_[col];
-          addToMatrix(row, col, -coefdelta);
+          addToMatrix(row, col, static_cast<double>(-coefdelta));
         }
       }
     };