Add Sturm sequence calculation

CMakeLists.txt

@@ -545,6 +545,7 @@ target_link_libraries(
   pugs
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsUtils
   PugsLanguage
   PugsLanguageAST
@@ -552,6 +553,7 @@ target_link_libraries(
   PugsLanguageAlgorithms
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsScheme
   PugsUtils
   PugsOutput

src/CMakeLists.txt

@@ -12,6 +12,9 @@ add_subdirectory(language)
 # Pugs algebra
 add_subdirectory(algebra)
 
+# Pugs analysis
+add_subdirectory(analysis)
+
 # Pugs mesh
 add_subdirectory(mesh)
 

src/analysis/CMakeLists.txt

+# ------------------- Source files --------------------
+
+add_library(
+  PugsAnalysis
+  SturmSequence.cpp
+)

src/analysis/SturmSequence.cpp

+#include <analysis/SturmSequence.hpp>
+
+SturmSequence::SturmSequence(const Polynomial1D& p0) : SturmSequence{Polynomial1D{p0}} {}
+
+SturmSequence::SturmSequence(Polynomial1D&& p0)
+{
+  m_sequence.reserve(p0.degree());
+  m_sequence.emplace_back(p0);
+  Polynomial1D p1 = derive(m_sequence[0]);
+  if (p1.degree() > 0 or (std::abs(p1.coefficient(0)) > 1E-14)) {
+    m_sequence.emplace_back(std::move(p1));
+
+    do {
+      Polynomial1D pk = -(m_sequence[m_sequence.size() - 2] % m_sequence[m_sequence.size() - 1]);
+      if (pk.degree() > 0 or (std::abs(pk.coefficient(0)) > 1E-14)) {
+        m_sequence.emplace_back(std::move(pk));
+      } else {
+        break;
+      }
+    } while (true);
+  }
+}
+
+size_t
+SturmSequence::numberOfSignChanges(const double x)
+{
+  std::vector<double> values;
+  values.reserve(m_sequence.size());
+  for (const auto& pk : m_sequence) {
+    values.emplace_back(pk(x));
+  }
+
+  size_t number_of_sign_changes = 0;
+  for (size_t i = 1; i < m_sequence.size(); ++i) {
+    number_of_sign_changes += (values[i - 1] * values[i] < 0);
+  }
+
+  return number_of_sign_changes;
+}

src/analysis/SturmSequence.hpp

+#ifndef STURM_SEQUENCE_HPP
+#define STURM_SEQUENCE_HPP
+
+#include <analysis/Polynomial1D.hpp>
+
+#include <vector>
+
+class SturmSequence
+{
+ private:
+  std::vector<Polynomial1D> m_sequence;
+
+ public:
+  PUGS_INLINE
+  size_t
+  size() const
+  {
+    return m_sequence.size();
+  }
+
+  size_t numberOfSignChanges(const double x);
+
+  SturmSequence(const Polynomial1D& p0);
+  SturmSequence(Polynomial1D&& p0);
+
+  SturmSequence(SturmSequence&&)      = delete;
+  SturmSequence(const SturmSequence&) = delete;
+
+  ~SturmSequence() = default;
+};
+
+#endif   // STURM_SEQUENCE_HPP

tests/CMakeLists.txt

@@ -89,6 +89,7 @@ add_executable (unit_tests
   test_PugsUtils.cpp
   test_RevisionInfo.cpp
   test_SparseMatrixDescriptor.cpp
+  test_SturmSequence.cpp
   test_SymbolTable.cpp
   test_Table.cpp
   test_Timer.cpp
@@ -128,6 +129,7 @@ target_link_libraries (unit_tests
   PugsLanguage
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsScheme
   PugsOutput
   PugsUtils
@@ -150,6 +152,7 @@ target_link_libraries (mpi_unit_tests
   PugsLanguageAlgorithms
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsUtils
   PugsLanguageUtils
   PugsScheme

tests/test_SturmSequence.cpp

+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/SturmSequence.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SturmSequence", "[analysis]")
+{
+  Polynomial1D p({-1, -1, 0, 1, 1});
+
+  SturmSequence sturm_sequence(p);
+
+  REQUIRE(sturm_sequence.size() == 5);
+
+  REQUIRE(sturm_sequence.numberOfSignChanges(-10) == 3);
+  REQUIRE(sturm_sequence.numberOfSignChanges(10) == 1);
+}