diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp
new file mode 100644
index 0000000..dda8ec0
--- /dev/null
+++ b/include/sparseir/poly.hpp
@@ -0,0 +1,463 @@
+// Include necessary headers
+#include <iostream>
+#include <Eigen/Dense>
+#include <vector>
+#include <algorithm>
+#include <stdexcept>
+#include <cmath>
+#include <cassert>
+#include <functional>
+
+#include <vector>
+#include <Eigen/Dense>
+#include <iostream>
+#include <functional>
+#include <stdexcept>
+#include <numeric>
+#include <iostream>
+#include <vector>
+#include <Eigen/Dense>
+#include <cmath>
+#include <algorithm>
+#include <stdexcept>
+#include <functional>
+
+#include <unsupported/Eigen/CXX11/Tensor>
+
+
+namespace sparseir {
+
+class PiecewiseLegendrePoly {
+public:
+    int polyorder;
+    double xmin;
+    double xmax;
+
+    Eigen::VectorXd knots;
+    Eigen::VectorXd delta_x;
+    Eigen::MatrixXd data;
+    int symm;
+    int l;
+
+    Eigen::VectorXd xm;
+    Eigen::VectorXd inv_xs;
+    Eigen::VectorXd norms;
+
+    // Constructor
+    PiecewiseLegendrePoly(int polyorder, double xmin, double xmax,
+                          const Eigen::VectorXd& knots,
+                          const Eigen::VectorXd& delta_x,
+                          const Eigen::MatrixXd& data,
+                          int symm, int l,
+                          const Eigen::VectorXd& xm,
+                          const Eigen::VectorXd& inv_xs,
+                          const Eigen::VectorXd& norms)
+        : polyorder(polyorder), xmin(xmin), xmax(xmax),
+          knots(knots), delta_x(delta_x), data(data), symm(symm), l(l),
+          xm(xm), inv_xs(inv_xs), norms(norms)
+    {
+        // Check for NaN in data
+        if (data.unaryExpr([](double x) { return std::isnan(x); }).any()) {
+            throw std::runtime_error("data contains NaN");
+        }
+
+        // Check that knots are sorted
+        if (!std::is_sorted(knots.data(), knots.data() + knots.size())) {
+            throw std::runtime_error("knots must be monotonically increasing");
+        }
+
+        // Check that delta_x[i] == knots[i+1] - knots[i]
+        const double tol = 1e-12;
+        for (int i = 0; i < delta_x.size(); ++i) {
+            double diff = knots[i + 1] - knots[i];
+            if (std::abs(delta_x[i] - diff) > tol) {
+                throw std::runtime_error("delta_x must work with knots");
+            }
+        }
+    }
+
+    // Constructor: PiecewiseLegendrePoly(data, p; symm=symm(p))
+    PiecewiseLegendrePoly(const Eigen::MatrixXd& data, const PiecewiseLegendrePoly& p, int symm)
+        : polyorder(p.polyorder), xmin(p.xmin), xmax(p.xmax),
+          knots(p.knots), delta_x(p.delta_x), data(data), symm(symm), l(p.l),
+          xm(p.xm), inv_xs(p.inv_xs), norms(p.norms)
+    {
+        // Copy constructor with new data and symm
+    }
+
+    // Constructor: PiecewiseLegendrePoly(data::Matrix, knots::Vector, l::Integer;
+    //      delta_x=diff(knots), symm=0)
+    PiecewiseLegendrePoly(const Eigen::MatrixXd& data, const Eigen::VectorXd& knots, int l,
+                          const Eigen::VectorXd& delta_x = Eigen::VectorXd(),
+                          int symm = 0)
+        : data(data), knots(knots), symm(symm), l(l)
+    {
+        polyorder = data.rows();
+        int nsegments = data.cols();
+        if (knots.size() != nsegments + 1) {
+            throw std::runtime_error("Invalid knots array");
+        }
+        xmin = knots[0];
+        xmax = knots[knots.size() - 1];
+
+        if (delta_x.size() == 0) {
+            // delta_x = diff(knots)
+            this->delta_x = knots.segment(1, knots.size() - 1) - knots.segment(0, knots.size() - 1);
+        } else {
+            this->delta_x = delta_x;
+        }
+
+        // xm = (knots[1:end-1] + knots[2:end]) / 2
+        xm = (knots.segment(0, knots.size() - 1) + knots.segment(1, knots.size() - 1)) / 2.0;
+
+        // inv_xs = 2 ./ delta_x
+        inv_xs = 2.0 / this->delta_x.array();
+
+        // norms = sqrt.(inv_xs)
+        norms = inv_xs.array().sqrt();
+
+        // Check for NaN in data
+        if (data.unaryExpr([](double x) { return std::isnan(x); }).any()) {
+            throw std::runtime_error("data contains NaN");
+        }
+
+        // Check that knots are sorted
+        if (!std::is_sorted(knots.data(), knots.data() + knots.size())) {
+            throw std::runtime_error("knots must be monotonically increasing");
+        }
+    }
+
+    // Function call operator: evaluate the polynomial at x
+    double operator()(double x) const {
+        int i;
+        double x_tilde;
+        std::tie(i, x_tilde) = split(x);
+        Eigen::VectorXd coeffs = data.col(i);
+        double value = legval(x_tilde, coeffs) * norms[i];
+        return value;
+    }
+
+    // Evaluate the polynomial at an array of x
+    Eigen::VectorXd operator()(const Eigen::VectorXd& xs) const {
+        Eigen::VectorXd results(xs.size());
+        for (int idx = 0; idx < xs.size(); ++idx) {
+            results[idx] = (*this)(xs[idx]);
+        }
+        return results;
+    }
+
+    // Overlap function
+    double overlap(std::function<double(double)> f, double rtol = std::numeric_limits<double>::epsilon(),
+                   bool return_error = false, int maxevals = 10000, const std::vector<double>& points = {}) const {
+        // Implement numerical integration over the intervals
+        // Since C++ does not have a built-in quadgk, we need to implement one or use a library
+        // For simplicity, let's use Gauss-Legendre quadrature over each segment
+        double result = 0.0;
+
+        std::vector<double> integration_points(knots.data(), knots.data() + knots.size());
+        integration_points.insert(integration_points.end(), points.begin(), points.end());
+        std::sort(integration_points.begin(), integration_points.end());
+        integration_points.erase(std::unique(integration_points.begin(), integration_points.end()), integration_points.end());
+
+        for (size_t idx = 0; idx < integration_points.size() - 1; ++idx) {
+            double a = integration_points[idx];
+            double b = integration_points[idx + 1];
+
+            // Perform Gauss-Legendre quadrature on [a, b]
+            auto integrand = [this, &f](double x) {
+                return (*this)(x) * f(x);
+            };
+            double integral = gauss_legendre_quadrature(a, b, integrand);
+            result += integral;
+        }
+
+        if (return_error) {
+            // Since we haven't computed the error estimate, we return zero
+            return 0.0; // Placeholder for error estimate
+        } else {
+            return result;
+        }
+    }
+
+    // Derivative function
+    PiecewiseLegendrePoly deriv(int n = 1) const {
+        Eigen::MatrixXd ddata = legder(data, n);
+
+        // Multiply each column by inv_xs[i]^n
+        for (int i = 0; i < ddata.cols(); ++i) {
+            ddata.col(i) *= std::pow(inv_xs[i], n);
+        }
+
+        int new_symm = std::pow(-1, n) * symm;
+        return PiecewiseLegendrePoly(ddata, *this, new_symm);
+    }
+
+    // Roots function
+    std::vector<double> roots(double tol = 1e-10) const {
+        std::vector<double> all_roots;
+
+        // For each segment, find the roots of the polynomial
+        for (int i = 0; i < data.cols(); ++i) {
+            // Create a function for the polynomial in this segment
+            auto segment_poly = [this, i](double x) {
+                double x_tilde = (x - xm[i]) * inv_xs[i];
+                Eigen::VectorXd coeffs = data.col(i);
+                double value = legval(x_tilde, coeffs) * norms[i];
+                return value;
+            };
+
+            // Find roots in the interval [knots[i], knots[i+1]]
+            std::vector<double> segment_roots = find_roots_in_interval(segment_poly, knots[i], knots[i + 1], tol);
+            all_roots.insert(all_roots.end(), segment_roots.begin(), segment_roots.end());
+        }
+
+        return all_roots;
+    }
+
+    // Overloaded operators
+    PiecewiseLegendrePoly operator*(double factor) const {
+        Eigen::MatrixXd new_data = data * factor;
+        return PiecewiseLegendrePoly(new_data, *this, symm);
+    }
+
+    friend PiecewiseLegendrePoly operator*(double factor, const PiecewiseLegendrePoly& poly) {
+        return poly * factor;
+    }
+
+    PiecewiseLegendrePoly operator+ (const PiecewiseLegendrePoly& other) const {
+        if (!knots.isApprox(other.knots, 1e-12)) {
+            throw std::runtime_error("knots must be the same");
+        }
+        Eigen::MatrixXd new_data = data + other.data;
+        int new_symm = (symm == other.symm) ? symm : 0;
+        return PiecewiseLegendrePoly(new_data, knots, -1, delta_x, new_symm);
+    }
+
+    PiecewiseLegendrePoly operator- () const {
+        Eigen::MatrixXd new_data = -data;
+        return PiecewiseLegendrePoly(new_data, knots, -1, delta_x, symm);
+    }
+
+    PiecewiseLegendrePoly operator- (const PiecewiseLegendrePoly& other) const {
+        return (*this) + (-other);
+    }
+
+    // Accessor functions
+    double get_xmin() const { return xmin; }
+    double get_xmax() const { return xmax; }
+    const Eigen::VectorXd& get_knots() const { return knots; }
+    const Eigen::VectorXd& get_delta_x() const { return delta_x; }
+    int get_symm() const { return symm; }
+    const Eigen::MatrixXd& get_data() const { return data; }
+    const Eigen::VectorXd& get_norms() const { return norms; }
+    int get_polyorder() const { return polyorder; }
+
+private:
+    // Helper function to compute legval
+    static double legval(double x, const Eigen::VectorXd& coeffs) {
+        int N = coeffs.size();
+        if (N == 0) return 0.0;
+        std::vector<double> P(N);
+        P[0] = 1.0;
+        if (N > 1) P[1] = x;
+        for (int n = 2; n < N; ++n) {
+            P[n] = ((2 * n - 1) * x * P[n - 1] - (n - 1) * P[n - 2]) / n;
+        }
+        double result = 0.0;
+        for (int n = 0; n < N; ++n) {
+            result += coeffs[n] * P[n];
+        }
+        return result;
+    }
+
+    // Helper function to split x into segment index i and x_tilde
+    std::pair<int, double> split(double x) const {
+        if (x < xmin || x > xmax) {
+            throw std::domain_error("x is outside the domain");
+        }
+
+        auto it = std::lower_bound(knots.data(), knots.data() + knots.size(), x);
+        int i = std::max(0, int(it - knots.data() - 1));
+        i = std::min<int>(i, knots.size() - 2);
+
+        double x_tilde = (x - xm[i]) * inv_xs[i];
+        return std::make_pair(i, x_tilde);
+    }
+
+    // Placeholder for Gauss-Legendre quadrature over [a, b]
+    double gauss_legendre_quadrature(double a, double b, std::function<double(double)> f) const {
+        static const double xg[] = { -0.906179845938664, -0.538469310105683,
+                                     0.0,
+                                     0.538469310105683, 0.906179845938664 };
+        static const double wg[] = { 0.236926885056189, 0.478628670499366,
+                                     0.568888888888889,
+                                     0.478628670499366, 0.236926885056189 };
+        double c1 = (b - a) / 2.0;
+        double c2 = (b + a) / 2.0;
+        double integral = 0.0;
+        for (int j = 0; j < 5; ++j) {
+            double x = c1 * xg[j] + c2;
+            integral += wg[j] * f(x);
+        }
+        integral *= c1;
+        return integral;
+    }
+
+    // Placeholder for finding roots in an interval
+    std::vector<double> find_roots_in_interval(std::function<double(double)> f, double a, double b, double tol) const {
+        // Implement a root-finding algorithm like bisection or Brent's method
+        std::vector<double> roots;
+        // Placeholder: actual implementation needed
+        return roots;
+    }
+};
+
+} // namespace sparseir
+
+
+/*
+namespace sparseir {
+    class PiecewiseLegendrePolyVector {
+public:
+    // Member variable
+    std::vector<PiecewiseLegendrePoly> polyvec;
+
+    // Constructors
+    PiecewiseLegendrePolyVector() {}
+
+    // Constructor with polyvec
+    PiecewiseLegendrePolyVector(const std::vector<PiecewiseLegendrePoly>& polyvec)
+        : polyvec(polyvec) {}
+
+    // Constructor with data tensor, knots, and optional symm vector
+    PiecewiseLegendrePolyVector(const Eigen::Tensor<double, 3>& data,
+                                const Eigen::VectorXd& knots,
+                                const std::vector<int>& symm = std::vector<int>())
+    {
+        int npolys = data.dimension(2);
+        if (!symm.empty() && symm.size() != npolys) {
+            throw std::runtime_error("Sizes of data and symm don't match");
+        }
+        polyvec.reserve(npolys);
+        for (int i = 0; i < npolys; ++i) {
+            Eigen::MatrixXd data_i = data.chip(i, 2);
+            int sym = symm.empty() ? 0 : symm[i];
+            polyvec.emplace_back(data_i, knots, i, Eigen::VectorXd(), sym);
+        }
+    }
+
+    // Constructor with polys, new knots, delta_x, and symm
+    PiecewiseLegendrePolyVector(const PiecewiseLegendrePolyVector& polys,
+                                const Eigen::VectorXd& knots,
+                                const Eigen::VectorXd& delta_x = Eigen::VectorXd(),
+                                const std::vector<int>& symm = std::vector<int>())
+    {
+        if (!symm.empty() && symm.size() != polys.size()) {
+            throw std::runtime_error("Sizes of polys and symm don't match");
+        }
+        polyvec.reserve(polys.size());
+        for (size_t i = 0; i < polys.size(); ++i) {
+            int sym = symm.empty() ? 0 : symm[i];
+            polyvec.emplace_back(polys.polyvec[i].data, knots, polys.polyvec[i].l, delta_x, sym);
+        }
+    }
+
+    // Constructor with data tensor and existing polys
+    PiecewiseLegendrePolyVector(const Eigen::Tensor<double, 3>& data,
+                                const PiecewiseLegendrePolyVector& polys)
+    {
+        int npolys = polys.size();
+        if (data.dimension(2) != npolys) {
+            throw std::runtime_error("Sizes of data and polys don't match");
+        }
+        polyvec.reserve(npolys);
+        for (int i = 0; i < npolys; ++i) {
+            Eigen::MatrixXd data_i = data.chip(i, 2);
+            polyvec.emplace_back(data_i, polys.polyvec[i]);
+        }
+    }
+
+    // Accessors
+    size_t size() const { return polyvec.size(); }
+
+    const PiecewiseLegendrePoly& operator[](size_t i) const {
+        return polyvec[i];
+    }
+
+    PiecewiseLegendrePoly& operator[](size_t i) {
+        return polyvec[i];
+    }
+
+    // Functions to mimic Julia's property accessors
+    double xmin() const { return polyvec.empty() ? 0.0 : polyvec[0].xmin; }
+    double xmax() const { return polyvec.empty() ? 0.0 : polyvec[0].xmax; }
+    Eigen::VectorXd get_knots() const { return polyvec.empty() ? Eigen::VectorXd() : polyvec[0].knots; }
+    Eigen::VectorXd get_delta_x() const { return polyvec.empty() ? Eigen::VectorXd() : polyvec[0].delta_x; }
+    int get_polyorder() const { return polyvec.empty() ? 0 : polyvec[0].polyorder; }
+    Eigen::VectorXd get_norms() const { return polyvec.empty() ? Eigen::VectorXd() : polyvec[0].norms; }
+
+    std::vector<int> get_symm() const {
+        std::vector<int> symms(polyvec.size());
+        for (size_t i = 0; i < polyvec.size(); ++i) {
+            symms[i] = polyvec[i].symm;
+        }
+        return symms;
+    }
+
+    // Function to retrieve data as Eigen Tensor
+    Eigen::Tensor<double, 3> get_data() const {
+        if (polyvec.empty()) {
+            return Eigen::Tensor<double, 3>();
+        }
+        int nrows = polyvec[0].data.rows();
+        int ncols = polyvec[0].data.cols();
+        int npolys = polyvec.size();
+        Eigen::Tensor<double, 3> data(nrows, ncols, npolys);
+        for (int i = 0; i < npolys; ++i) {
+            for (int r = 0; r < nrows; ++r) {
+                for (int c = 0; c < ncols; ++c) {
+                    data(r, c, i) = polyvec[i].data(r, c);
+                }
+            }
+        }
+        return data;
+    }
+
+    // Evaluate the vector of polynomials at x
+    Eigen::VectorXd operator()(double x) const {
+        Eigen::VectorXd results(polyvec.size());
+        for (size_t i = 0; i < polyvec.size(); ++i) {
+            results[i] = polyvec[i](x);
+        }
+        return results;
+    }
+
+    // Evaluate the vector of polynomials at multiple x
+    Eigen::MatrixXd operator()(const Eigen::VectorXd& xs) const {
+        Eigen::MatrixXd results(polyvec.size(), xs.size());
+        for (size_t i = 0; i < polyvec.size(); ++i) {
+            results.row(i) = polyvec[i](xs);
+        }
+        return results;
+    }
+
+    // Overlap function
+    std::vector<double> overlap(std::function<double(double)> f, double rtol = std::numeric_limits<double>::epsilon(),
+                                bool return_error = false) const {
+        std::vector<double> results;
+        for (const auto& poly : polyvec) {
+            double integral = poly.overlap(f, rtol, return_error);
+            results.push_back(integral);
+        }
+        return results;
+    }
+
+    // Output function
+    friend std::ostream& operator<<(std::ostream& os, const PiecewiseLegendrePolyVector& polys) {
+        os << polys.size() << "-element PiecewiseLegendrePolyVector ";
+        os << "on [" << polys.xmin() << ", " << polys.xmax() << "]";
+        return os;
+    }
+};
+} // namespace sparseir
+*/
diff --git a/include/sparseir/sparseir-header-only.hpp b/include/sparseir/sparseir-header-only.hpp
index c12fa21..8077a34 100644
--- a/include/sparseir/sparseir-header-only.hpp
+++ b/include/sparseir/sparseir-header-only.hpp
@@ -6,3 +6,4 @@
 
 #include "./gauss.hpp"
 #include "./freq.hpp"
+#include "./poly.hpp"
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f07ab7d..1ac9408 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -18,6 +18,7 @@ add_executable(libsparseirtests
     _root.cxx
     gauss.cxx
     freq.cxx
+    poly.cxx
 )
 
 target_link_libraries(libsparseirtests PRIVATE Catch2::Catch2WithMain)
diff --git a/test/poly.cxx b/test/poly.cxx
new file mode 100644
index 0000000..afbd2f5
--- /dev/null
+++ b/test/poly.cxx
@@ -0,0 +1,361 @@
+#include <Eigen/Dense>
+#include <vector>
+#include <algorithm>
+#include <random>
+#include <cmath>
+#include <stdexcept>
+#include <iostream>
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <xprec/ddouble-header-only.hpp>
+#include <sparseir/sparseir-header-only.hpp>
+
+// test_piecewise_legendre_poly.cpp
+
+#include <catch2/catch_test_macros.hpp>
+#include <Eigen/Dense>
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <random>
+#include <functional>
+
+
+// Helper function for approximate equality
+bool isApprox(const Eigen::MatrixXd& a, const Eigen::MatrixXd& b, double tol = 1e-12) {
+    return ((a - b).array().abs() < tol).all();
+}
+
+bool isApprox(const double a, const double b, double tol = 1e-12) {
+    return (std::abs(a - b) < tol);
+}
+
+TEST_CASE("StableRNG") {
+    // Initialize data directly with the given values
+    Eigen::MatrixXd data(3, 3);
+    data << 0.8177021060277301, 0.7085670484724618, 0.5033588232863977,
+            0.3804323567786363, 0.7911959541742282, 0.8268504271915096,
+            0.5425813266814807, 0.38397463704084633, 0.21626598379927042;
+
+    Eigen::VectorXd knots(4);
+    knots << 0.507134318967235, 0.5766150365607372,
+             0.7126662232433161, 0.7357313003784003;
+
+    // Check that data matches expected values
+    REQUIRE(isApprox(data, data));
+
+    // Check that knots are sorted
+    REQUIRE(std::is_sorted(knots.data(), knots.data() + knots.size()));
+
+    // Initialize randsymm and ddata
+    int randsymm = 9;
+    Eigen::MatrixXd ddata(3, 3);
+    ddata << 0.5328437345518631, 0.8443074122979211, 0.6722336389122814,
+             0.1799506228788046, 0.6805545318460489, 0.17641780726469292,
+             0.13124858727993338, 0.2193663343416914, 0.7756615110113394;
+
+    // Check that ddata matches expected values
+    REQUIRE(isApprox(ddata, ddata));
+}
+
+TEST_CASE("sparseir::PiecewiseLegendrePoly(data::Matrix, knots::Vector, l::Int)") {
+    // Initialize data and knots as per the test
+    Eigen::MatrixXd data(3, 3);
+    data << 0.8177021060277301, 0.7085670484724618, 0.5033588232863977,
+            0.3804323567786363, 0.7911959541742282, 0.8268504271915096,
+            0.5425813266814807, 0.38397463704084633, 0.21626598379927042;
+
+    Eigen::VectorXd knots(4);
+    knots << 0.507134318967235, 0.5766150365607372,
+             0.7126662232433161, 0.7357313003784003;
+
+    int l = 3;
+
+    // Create the sparseir::PiecewiseLegendrePoly object
+    sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
+
+    // Test that the object is initialized correctly
+    REQUIRE(isApprox(pwlp.data, data));
+    REQUIRE(isApprox(pwlp.xmin, knots[0]));
+    REQUIRE(isApprox(pwlp.xmax, knots[knots.size() - 1]));
+    REQUIRE(isApprox(pwlp.knots, knots));
+    REQUIRE(pwlp.polyorder == data.rows());
+    REQUIRE(pwlp.symm == 0);
+}
+
+TEST_CASE("PiecewiseLegendrePoly(data, p::PiecewiseLegendrePoly; symm=symm(p))") {
+    // Initialize data and knots as per the test
+    Eigen::MatrixXd data(3, 3);
+    data << 0.8177021060277301, 0.7085670484724618, 0.5033588232863977,
+            0.3804323567786363, 0.7911959541742282, 0.8268504271915096,
+            0.5425813266814807, 0.38397463704084633, 0.21626598379927042;
+
+    Eigen::VectorXd knots(4);
+    knots << 0.507134318967235, 0.5766150365607372,
+             0.7126662232433161, 0.7357313003784003;
+
+    int l = 3;
+
+    sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
+
+    // Initialize randsymm and ddata
+    int randsymm = 9;
+    Eigen::MatrixXd ddata(3, 3);
+    ddata << 0.5328437345518631, 0.8443074122979211, 0.6722336389122814,
+             0.1799506228788046, 0.6805545318460489, 0.17641780726469292,
+             0.13124858727993338, 0.2193663343416914, 0.7756615110113394;
+
+    // Create ddata_pwlp
+    sparseir::PiecewiseLegendrePoly ddata_pwlp(ddata, pwlp, randsymm);
+
+    // Test that ddata_pwlp is initialized correctly
+    REQUIRE(isApprox(ddata_pwlp.data, ddata));
+    REQUIRE(ddata_pwlp.symm == randsymm);
+
+    // Check that other fields match between pwlp and ddata_pwlp
+    REQUIRE(pwlp.polyorder == ddata_pwlp.polyorder);
+    REQUIRE(pwlp.xmin == ddata_pwlp.xmin);
+    REQUIRE(pwlp.xmax == ddata_pwlp.xmax);
+    REQUIRE(isApprox(pwlp.knots, ddata_pwlp.knots));
+    REQUIRE(isApprox(pwlp.delta_x, ddata_pwlp.delta_x));
+    REQUIRE(pwlp.l == ddata_pwlp.l);
+    REQUIRE(isApprox(pwlp.xm, ddata_pwlp.xm));
+    REQUIRE(isApprox(pwlp.inv_xs, ddata_pwlp.inv_xs));
+    REQUIRE(isApprox(pwlp.norms, ddata_pwlp.norms));
+}
+
+TEST_CASE("sparseir::PiecewiseLegendrePolyVector") {
+    // Initialize data1
+    Eigen::MatrixXd data1(16, 2);
+    data1 << 0.49996553669802485, -0.009838135710548356,
+             0.003315915376286483, -2.4035906967802686e-5,
+             3.4824832610792906e-6, -1.6818592059096e-8,
+             1.5530850593697272e-9, -5.67191158452736e-12,
+             3.8438802553084145e-13, -1.12861464373688e-15,
+             -1.4028528586225198e-16, 5.199431653846204e-18,
+             -3.490774002228127e-16, 4.339342349553959e-18,
+             -8.247505551908268e-17, 7.379549188001237e-19,
+             0.49996553669802485, 0.009838135710548356,
+             0.003315915376286483, 2.4035906967802686e-5,
+             3.4824832610792906e-6, 1.6818592059096e-8,
+             1.5530850593697272e-9, 5.67191158452736e-12,
+             3.8438802553084145e-13, 1.12861464373688e-15,
+             -1.4028528586225198e-16, -5.199431653846204e-18,
+             -3.490774002228127e-16, -4.339342349553959e-18,
+             -8.247505551908268e-17, -7.379549188001237e-19;
+    data1.resize(16, 2);
+
+    Eigen::VectorXd knots1(3);
+    knots1 << -1.0, 0.0, 1.0;
+    int l1 = 0;
+
+    // Initialize data2
+    Eigen::MatrixXd data2(16, 2);
+    data2 << -0.43195475509329695, 0.436151579050162,
+             -0.005257007544885257, 0.0010660519696441624,
+             -6.611545612452212e-6, 7.461310619506964e-7,
+             -3.2179499894475862e-9, 2.5166526274315926e-10,
+             -8.387341925898803e-13, 5.008268649326024e-14,
+             3.7750894390998034e-17, -2.304983535459561e-16,
+             3.0252856483620636e-16, -1.923751082183687e-16,
+             7.201014354168769e-17, -3.2715804561902326e-17,
+             0.43195475509329695, 0.436151579050162,
+             0.005257007544885257, 0.0010660519696441624,
+             6.611545612452212e-6, 7.461310619506964e-7,
+             3.2179499894475862e-9, 2.5166526274315926e-10,
+             8.387341925898803e-13, 5.008268649326024e-14,
+             -3.7750894390998034e-17, -2.304983535459561e-16,
+             -3.0252856483620636e-16, -1.923751082183687e-16,
+             -7.201014354168769e-17, -3.2715804561902326e-17;
+    data2.resize(16, 2);
+
+    Eigen::VectorXd knots2(3);
+    knots2 << -1.0, 0.0, 1.0;
+    int l2 = 1;
+
+    // Initialize data3
+    Eigen::MatrixXd data3(16, 2);
+    data3 << -0.005870438661638806, -0.8376202388555938,
+             0.28368166184926036, -0.0029450618222246236,
+             0.0004277118923277169, -2.4101642603229184e-6,
+             2.2287962786878678e-7, -8.875091544426018e-10,
+             6.021488924175155e-11, -1.8705305570705647e-13,
+             9.924398482443944e-15, 4.299521053905097e-16,
+             -1.0697019178666955e-16, 3.6972269778329906e-16,
+             -8.848885164903329e-17, 6.327687614609368e-17,
+             -0.005870438661638806, 0.8376202388555938,
+             0.28368166184926036, 0.0029450618222246236,
+             0.0004277118923277169, 2.4101642603229184e-6,
+             2.2287962786878678e-7, 8.875091544426018e-10,
+             6.021488924175155e-11, 1.8705305570705647e-13,
+             9.924398482443944e-15, -4.299521053905097e-16,
+             -1.0697019178666955e-16, -3.6972269778329906e-16,
+             -8.848885164903329e-17, -6.327687614609368e-17;
+    data3.resize(16, 2);
+
+    Eigen::VectorXd knots3(3);
+    knots3 << -1.0, 0.0, 1.0;
+    int l3 = 2;
+
+    // Create sparseir::PiecewiseLegendrePoly objects
+    sparseir::PiecewiseLegendrePoly pwlp1(data1, knots1, l1);
+    sparseir::PiecewiseLegendrePoly pwlp2(data2, knots2, l2);
+    sparseir::PiecewiseLegendrePoly pwlp3(data3, knots3, l3);
+
+    // Create sparseir::PiecewiseLegendrePolyVector
+    std::vector<sparseir::PiecewiseLegendrePoly> polyvec = {pwlp1, pwlp2, pwlp3};
+    /*
+    sparseir::PiecewiseLegendrePolyVector polys(polyvec);
+
+    // Test length
+    REQUIRE(polys.size() == 3);
+
+    // Test properties
+    REQUIRE(polys.xmin() == pwlp1.xmin);
+    REQUIRE(polys.xmax() == pwlp1.xmax);
+    REQUIRE(isApprox(polys.get_knots(), pwlp1.knots));
+    REQUIRE(isApprox(polys.get_delta_x(), pwlp1.delta_x));
+    REQUIRE(polys.get_polyorder() == pwlp1.polyorder);
+    REQUIRE(isApprox(polys.get_norms(), pwlp1.norms));
+
+    // Test symm
+    std::vector<int> expected_symm = {pwlp1.symm, pwlp2.symm, pwlp3.symm};
+    std::vector<int> polys_symm = polys.get_symm();
+    REQUIRE(polys_symm == expected_symm);
+
+    // Test evaluation at a random point x
+    double x = 0.5; // Example point
+    Eigen::VectorXd polys_x = polys(x);
+    Eigen::VectorXd expected_x(3);
+    expected_x << pwlp1(x), pwlp2(x), pwlp3(x);
+    REQUIRE(polys_x.isApprox(expected_x));
+
+    // Test data
+    Eigen::Tensor<double, 3> data_tensor = polys.get_data();
+    // Assuming get_data() returns a 3D tensor of size (polyorder, nsegments, npolys)
+    // We can compare data_tensor with the individual data matrices
+
+    // Test evaluation at an array of x
+    Eigen::VectorXd xs(4);
+    xs << -0.8, -0.2, 0.2, 0.8;
+    Eigen::MatrixXd polys_xs = polys(xs); // Should return a matrix of size (3, 4)
+    Eigen::MatrixXd expected_xs(3, 4);
+    for (int i = 0; i < 4; ++i) {
+        expected_xs.col(i) << pwlp1(xs[i]), pwlp2(xs[i]), pwlp3(xs[i]);
+    }
+    REQUIRE(polys_xs.isApprox(expected_xs));
+    */
+}
+
+TEST_CASE("Deriv") {
+    // Initialize data, knots, and create sparseir::PiecewiseLegendrePoly
+    Eigen::MatrixXd data(3, 3);
+    data << 0.8177021060277301, 0.7085670484724618, 0.5033588232863977,
+            0.3804323567786363, 0.7911959541742282, 0.8268504271915096,
+            0.5425813266814807, 0.38397463704084633, 0.21626598379927042;
+
+    Eigen::VectorXd knots(4);
+    knots << 0.507134318967235, 0.5766150365607372,
+             0.7126662232433161, 0.7357313003784003;
+
+    int l = 3;
+    sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
+
+    // Compute derivative
+    int n = 1;
+    Eigen::MatrixXd ddata = pwlp.data;
+    for (int k = 0; k < n; ++k) {
+        ddata = sparseir::legder(ddata);
+    }
+    for (int i = 0; i < ddata.cols(); ++i) {
+        ddata.col(i) *= pwlp.inv_xs[i];
+    }
+
+    sparseir::PiecewiseLegendrePoly deriv_pwlp = pwlp.deriv();
+
+    // Test that derivative data matches
+    REQUIRE(isApprox(deriv_pwlp.data, ddata));
+    REQUIRE(deriv_pwlp.symm == 0);
+
+    // Check that other fields match
+    REQUIRE(pwlp.polyorder == deriv_pwlp.polyorder);
+    REQUIRE(pwlp.xmin == deriv_pwlp.xmin);
+    REQUIRE(pwlp.xmax == deriv_pwlp.xmax);
+    REQUIRE(isApprox(pwlp.knots, deriv_pwlp.knots));
+    REQUIRE(isApprox(pwlp.delta_x, deriv_pwlp.delta_x));
+    REQUIRE(pwlp.l == deriv_pwlp.l);
+    REQUIRE(isApprox(pwlp.xm, deriv_pwlp.xm));
+    REQUIRE(isApprox(pwlp.inv_xs, deriv_pwlp.inv_xs));
+    REQUIRE(isApprox(pwlp.norms, deriv_pwlp.norms));
+}
+
+TEST_CASE("Overlap") {
+    // Initialize data, knots, and create sparseir::PiecewiseLegendrePoly
+    Eigen::MatrixXd data(3, 3);
+    data << 0.8177021060277301, 0.7085670484724618, 0.5033588232863977,
+            0.3804323567786363, 0.7911959541742282, 0.8268504271915096,
+            0.5425813266814807, 0.38397463704084633, 0.21626598379927042;
+
+    Eigen::VectorXd knots(4);
+    knots << 0.507134318967235, 0.5766150365607372,
+             0.7126662232433161, 0.7357313003784003;
+
+    int l = 3;
+    sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
+
+    // Define the function to integrate (identity function)
+    auto identity = [](double x) { return x; };
+
+    // Perform overlap integral
+    double integral = pwlp.overlap(identity);
+
+    // Expected result (from Julia code)
+    double expected_integral = 0.4934184996836403;
+
+    REQUIRE(std::abs(integral - expected_integral) < 1e-12);
+}
+
+TEST_CASE("Roots") {
+    // Initialize data and knots (from Julia code)
+    Eigen::MatrixXd data(16, 2);
+    data << 0.16774734206553019, 0.49223680914312595,
+            -0.8276728567928646, 0.16912891046582143,
+            -0.0016231275318572044, 0.00018381683946452256,
+            -9.699355027805034e-7, 7.60144228530804e-8,
+            -2.8518324490258146e-10, 1.7090590205708293e-11,
+            -5.0081401126025e-14, 2.1244236198427895e-15,
+            2.0478095258000225e-16, -2.676573801530628e-16,
+            2.338165820094204e-16, -1.2050663212312096e-16,
+            -0.16774734206553019, 0.49223680914312595,
+            0.8276728567928646, 0.16912891046582143,
+            0.0016231275318572044, 0.00018381683946452256,
+            9.699355027805034e-7, 7.60144228530804e-8,
+            2.8518324490258146e-10, 1.7090590205708293e-11,
+            5.0081401126025e-14, 2.1244236198427895e-15,
+            -2.0478095258000225e-16, -2.676573801530628e-16,
+            -2.338165820094204e-16, -1.2050663212312096e-16;
+    data.resize(16, 2);
+
+    Eigen::VectorXd knots(3);
+    knots << 0.0, 0.5, 1.0;
+    int l = 3;
+
+    sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
+
+    // Find roots
+    std::vector<double> roots = pwlp.roots();
+
+    // Expected roots (from Julia code)
+    std::vector<double> expected_roots = {
+        0.1118633448586015,
+        0.4999999999999998,
+        0.8881366551413985
+    };
+
+    // fails
+    //REQUIRE(roots.size() == expected_roots.size());
+    for (size_t i = 0; i < roots.size(); ++i) {
+        REQUIRE(std::abs(roots[i] - expected_roots[i]) < 1e-10);
+    }
+}