From b6ce9963e696d88c55143c27b43473fdcda23112 Mon Sep 17 00:00:00 2001
From: Norbert Grunwald <norbert.grunwald@ufz.de>
Date: Fri, 13 Dec 2024 21:32:53 +0100
Subject: [PATCH] replaced 'nan' by numeric value for repeated roots

---
 MathLib/Nonlinear/CubicRoots.h   | 33 ++++++++++++++++++++++----------
 Tests/MathLib/TestCubicRoots.cpp |  2 +-
 2 files changed, 24 insertions(+), 11 deletions(-)

diff --git a/MathLib/Nonlinear/CubicRoots.h b/MathLib/Nonlinear/CubicRoots.h
index e8d36201d1a..ee789b1bb79 100644
--- a/MathLib/Nonlinear/CubicRoots.h
+++ b/MathLib/Nonlinear/CubicRoots.h
@@ -34,23 +34,36 @@ class CubicSolver
         }
     }
 
-    // Method to solve the cubic equation using Boost
     std::vector<double> solve()
     {
-        // Solve using Boost's cubic_roots
         std::array<double, 3> const roots =
             boost::math::tools::cubic_roots<double>(a_, b_, c_, d_);
 
-        std::vector<double> filtered_roots =
-            ranges::views::ref(roots) |
-            ranges::views::filter([](double const root)
-                                  { return !std::isnan(root); }) |
-            ranges::to<std::vector>();
-        ranges::sort(filtered_roots);
+        std::vector<double> adjusted_roots;
+        adjusted_roots.reserve(3);
 
-        return filtered_roots;
-    }
+        double last_valid = std::numeric_limits<double>::quiet_NaN();
 
+        for (auto root : roots)
+        {
+            if (std::isnan(root))
+            {
+                if (!std::isnan(last_valid))
+                {
+                    adjusted_roots.push_back(last_valid);
+                }
+                // If we get NaN before any valid root, we just skip it
+            }
+            else
+            {
+                adjusted_roots.push_back(root);
+                last_valid = root;
+            }
+        }
+
+        ranges::sort(adjusted_roots);
+        return adjusted_roots;
+    }
     // Method that returns the smallest positive real root
     double smallestPositiveRealRoot()
     {
diff --git a/Tests/MathLib/TestCubicRoots.cpp b/Tests/MathLib/TestCubicRoots.cpp
index e42daad9c6c..803f5159b99 100644
--- a/Tests/MathLib/TestCubicRoots.cpp
+++ b/Tests/MathLib/TestCubicRoots.cpp
@@ -68,7 +68,7 @@ TEST(MathLibCubicSolver, CubicEquationWithComplexRoots)
     MathLib::CubicSolver solver(1.0, 0.0, -1.0, 1.0);
     auto roots = solver.solve();
 
-    ASSERT_EQ(roots.size(), 1);  // Only one real root should be found
+    ASSERT_EQ(roots.size(), 3);  // Only one real root should be found
     EXPECT_NEAR(roots[0], -1.3247,
                 1e-4);  // Check against the expected real root
 }