diff --git a/src/analysis/SturmSequence.cpp b/src/analysis/SturmSequence.cpp
index bda7f5fa58be9bfcee2c2ef1b207583053ba8eb5..f4af48ce852d5afb3f4f2942fb4bc1a66ac6da4f 100644
--- a/src/analysis/SturmSequence.cpp
+++ b/src/analysis/SturmSequence.cpp
@@ -22,7 +22,7 @@ SturmSequence::SturmSequence(Polynomial1D&& p0)
 }
 
 size_t
-SturmSequence::numberOfSignChanges(const double x)
+SturmSequence::numberOfSignChanges(const double x) const
 {
   std::vector<double> values;
   values.reserve(m_sequence.size());
diff --git a/src/analysis/SturmSequence.hpp b/src/analysis/SturmSequence.hpp
index 5445cb65f96887faca185f713e8cac43c93fc312..e7377ffcf586c6003650c8ffa76c9a25ed804533 100644
--- a/src/analysis/SturmSequence.hpp
+++ b/src/analysis/SturmSequence.hpp
@@ -11,6 +11,14 @@ class SturmSequence
   std::vector<Polynomial1D> m_sequence;
 
  public:
+  PUGS_INLINE
+  const Polynomial1D&
+  p(const size_t i) const
+  {
+    Assert(i < m_sequence.size());
+    return m_sequence[i];
+  }
+
   PUGS_INLINE
   size_t
   size() const
@@ -18,7 +26,7 @@ class SturmSequence
     return m_sequence.size();
   }
 
-  size_t numberOfSignChanges(const double x);
+  size_t numberOfSignChanges(const double x) const;
 
   SturmSequence(const Polynomial1D& p0);
   SturmSequence(Polynomial1D&& p0);