Fix

include/sparseir/_root.hpp

@@ -235,7 +235,7 @@ std::vector<double> discrete_extrema(F f, const std::vector<double> &xgrid)
     std::vector<double> res;
     for (size_t i = 0; i < a.size(); ++i) {
         // abs ∘ f
-        auto abf = [f](double x) { return std::abs(f(x)); };
+        auto abf = [f](double x) { return std::fabs(f(x)); };
         res.push_back(
             bisect_discr_extremum(abf, a[i], b[i], absf_a[i], absf_b[i]));
     }

include/sparseir/poly.hpp

@@ -250,8 +250,6 @@ class PiecewiseLegendrePoly {
     {
         Eigen::VectorXd grid = this->knots;
 
-        std::cout << "grid: " << grid.transpose() << "\n";
-
         Eigen::VectorXd refined_grid = refine_grid(grid, 2);
         auto f = [this](double x) { return this->operator()(x); };
         // convert to std::vector<double>

test/poly.cxx

@@ -329,10 +329,13 @@ TEST_CASE("Roots")
 
     REQUIRE(roots.size() == expected_roots.size());
     for(Eigen::Index i = 0; i < roots.size(); ++i) {
-        REQUIRE(std::abs(roots[i] - expected_roots[i]) < 1e-10);
-        //Verify roots are in domain
+        REQUIRE(std::fabs(roots[i] - expected_roots[i]) < 1e-10);
+        // Verifying the polynomial evaluates to zero
+        //at the roots with tight tolerance
+        REQUIRE(std::fabs(pwlp(roots[i])) < 1e-10);
+        // Verify roots are in domain
         REQUIRE(roots[i] >= knots[0]);
-        REQUIRE(roots[i] <= knots[knots.size()-1]);
+        REQUIRE(roots[i] <= knots[knots.size() - 1]);
         // Verify these are actually roots
         REQUIRE(std::abs(pwlp(roots[i])) < 1e-10);
     }