diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b3521136bb4442a17573ebd06a7d0c34d8cda685
--- /dev/null
+++ b/src/analysis/PolynomialP.hpp
@@ -0,0 +1,492 @@
+#ifndef POLYNOMIALP_HPP
+#define POLYNOMIALP_HPP
+
+#include <algebra/TinyVector.hpp>
+
+template <size_t N, size_t Dimension>
+class PolynomialP
+{
+ private:
+ static constexpr size_t size_coef = [] {
+ if constexpr (Dimension == 1) {
+ return N + 1;
+ } else if constexpr (Dimension == 2) {
+ return (N + 1) * (N + 2) / 2;
+ } else {
+ static_assert(Dimension == 3);
+ return (N + 1) * (N + 2) * (N + 3) / 6;
+ }
+ }();
+
+ using Coefficients = TinyVector<size_coef, double>;
+ Coefficients m_coefficients;
+ static_assert((N >= 0), "PolynomialP degree must be non-negative");
+ static_assert((Dimension > 0), "PolynomialP dimension must be positive");
+ static_assert((Dimension <= 3), "PolynomialP dimension must no greater than three");
+
+ public:
+ PUGS_INLINE
+ constexpr size_t
+ degree() const
+ {
+ return N;
+ }
+
+ constexpr size_t
+ dim() const
+ {
+ return Dimension;
+ }
+
+ PUGS_INLINE
+ constexpr const TinyVector<size_coef, double>&
+ coefficients() const
+ {
+ return m_coefficients;
+ }
+
+ PUGS_INLINE
+ constexpr TinyVector<size_coef, double>&
+ coefficients()
+ {
+ return m_coefficients;
+ }
+
+ PUGS_INLINE constexpr bool
+ operator==(const PolynomialP& q) const
+ {
+ return m_coefficients == q.m_coefficients;
+ }
+
+ PUGS_INLINE constexpr PolynomialP(const TinyVector<size_coef, double>& coefficients) noexcept
+ : m_coefficients{coefficients}
+ {}
+
+ PUGS_INLINE
+ constexpr PolynomialP(TinyVector<size_coef, double>&& coefficients) noexcept : m_coefficients{coefficients} {}
+
+ PUGS_INLINE
+ constexpr bool
+ operator!=(const PolynomialP& q) const
+ {
+ return not this->operator==(q);
+ }
+
+ PUGS_INLINE constexpr PolynomialP
+ operator+(const PolynomialP Q) const
+ {
+ PolynomialP<N, Dimension> P(m_coefficients);
+ for (size_t i = 0; i < size_coef; ++i) {
+ P.coefficients()[i] += Q.coefficients()[i];
+ }
+ return P;
+ }
+
+ PUGS_INLINE constexpr PolynomialP
+ operator-() const
+ {
+ PolynomialP<N, Dimension> P;
+ P.coefficients() = -coefficients();
+ return P;
+ }
+
+ PUGS_INLINE constexpr PolynomialP() noexcept = default;
+ ~PolynomialP() = default;
+};
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<std::max(M, N)>
+// operator-(const PolynomialP<M>& Q) const
+// {
+// PolynomialP<std::max(M, N)> P;
+// if constexpr (M > N) {
+// P.coefficients() = -Q.coefficients();
+// for (size_t i = 0; i <= N; ++i) {
+// P.coefficients()[i] += coefficients()[i];
+// }
+// } else {
+// P.coefficients() = coefficients();
+// for (size_t i = 0; i <= M; ++i) {
+// P.coefficients()[i] -= Q.coefficients()[i];
+// }
+// }
+// return P;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<N>&
+// operator=(const PolynomialP<M>& Q)
+// {
+// coefficients() = zero;
+// for (size_t i = N + 1; i <= M; ++i) {
+// Assert(Q.coefficients()[i] == 0, "degree of polynomialP to small in assignation");
+// }
+// // static_assert(N >= M, "degree of polynomialP to small in assignation");
+// for (size_t i = 0; i <= std::min(M, N); ++i) {
+// coefficients()[i] = Q.coefficients()[i];
+// }
+// return *this;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<N>&
+// operator+=(const PolynomialP<M>& Q)
+// {
+// static_assert(N >= M, "PolynomialP degree to small in affectation addition");
+// for (size_t i = 0; i <= M; ++i) {
+// coefficients()[i] += Q.coefficients()[i];
+// }
+// return *this;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<N>&
+// operator-=(const PolynomialP<M>& Q)
+// {
+// static_assert(N >= M, "PolynomialP degree to small in affectation addition");
+// for (size_t i = 0; i <= M; ++i) {
+// coefficients()[i] -= Q.coefficients()[i];
+// }
+// return *this;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<M + N>
+// operator*(const PolynomialP<M>& Q) const
+// {
+// PolynomialP<M + N> P;
+// P.coefficients() = zero;
+// for (size_t i = 0; i <= N; ++i) {
+// for (size_t j = 0; j <= M; ++j) {
+// P.coefficients()[i + j] += coefficients()[i] * Q.coefficients()[j];
+// }
+// }
+// return P;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<N>&
+// operator*=(const PolynomialP<M>& Q)
+// {
+// static_assert(N >= M, "Degree to small in affectation product");
+// for (size_t i = N - M + 1; i <= N; ++i) {
+// Assert(coefficients()[i] == 0, "Degree of affectation product greater than the degree of input polynomialP");
+// }
+// PolynomialP<N> P(zero);
+// for (size_t i = 0; i <= N - M; ++i) {
+// for (size_t j = 0; j <= M; ++j) {
+// P.coefficients()[i + j] += coefficients()[i] * Q.coefficients()[j];
+// }
+// }
+// coefficients() = P.coefficients();
+// return *this;
+// }
+
+// PUGS_INLINE
+// constexpr PolynomialP<N>
+// operator*(const double& lambda) const
+// {
+// TinyVector<N + 1> mult_coefs = lambda * m_coefficients;
+// PolynomialP<N> M(mult_coefs);
+// return M;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<M * N>
+// compose(const PolynomialP<M>& Q) const
+// {
+// PolynomialP<M * N> P;
+// P.coefficients() = zero;
+// PolynomialP<M * N> R;
+// R.coefficients() = zero;
+
+// for (size_t i = 0; i <= N; ++i) {
+// R = Q.template pow<N>(i) * coefficients()[i];
+// P += R; // R;
+// }
+// return P;
+// }
+
+// template <size_t M, size_t I>
+// PUGS_INLINE constexpr PolynomialP<M * N>
+// power(const PolynomialP<M>& Q) const
+// {
+// return coefficients()[I] * Q.template pow<N>(I);
+// }
+
+// template <size_t M, size_t... I>
+// PUGS_INLINE constexpr PolynomialP<M * N>
+// compose2(const PolynomialP<M>& Q, std::index_sequence<I...>) const
+// {
+// PolynomialP<M * N> P;
+// P.coefficients() = zero;
+// P = (power<M, I>(Q) + ...);
+// return P;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<M * N>
+// compose2(const PolynomialP<M>& Q) const
+// {
+// PolynomialP<M * N> P;
+// P.coefficients() = zero;
+// using IndexSequence = std::make_index_sequence<N + 1>;
+// return compose2<M>(Q, IndexSequence{});
+// // for (size_t i = 0; i <= N; ++i) {
+// // P += Q.template pow<N>(i) * coefficients()[i];
+// // }
+// return P;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<M * N>
+// operator()(const PolynomialP<M>& Q) const
+// {
+// PolynomialP<M * N> P;
+// P.coefficients() = zero;
+// PolynomialP<M * N> R;
+// R.coefficients() = zero;
+
+// for (size_t i = 0; i <= N; ++i) {
+// R = Q.template pow<N>(i) * coefficients()[i];
+// P += R; // R;
+// }
+// return P;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr PolynomialP<M * N>
+// pow(size_t power) const
+// {
+// Assert(power <= M, "You declared a polynomialP of degree too small for return of the pow function");
+// PolynomialP<M * N> R;
+// R.coefficients() = zero;
+// if (power == 0) {
+// R.coefficients()[0] = 1;
+// } else {
+// R = *this;
+// for (size_t i = 1; i < power; ++i) {
+// R = R * *this;
+// }
+// }
+// return R;
+// }
+
+// PUGS_INLINE
+// constexpr friend PolynomialP<N>
+// operator*(const double& lambda, const PolynomialP<N> P)
+// {
+// return P * lambda;
+// }
+// // evaluation using Horner's method https://en.wikipedia.org/wiki/Horner's_method
+// PUGS_INLINE
+// constexpr double
+// evaluate(const double& x) const
+// {
+// TinyVector<N + 1> coefs = this->coefficients();
+// double bcoef = coefs[N];
+// for (size_t i = N; i > 0; --i) {
+// bcoef *= x;
+// bcoef += coefs[i - 1];
+// }
+// return bcoef;
+// }
+
+// PUGS_INLINE
+// constexpr double
+// operator()(const double x) const
+// {
+// TinyVector<N + 1> coefs = this->coefficients();
+// double bcoef = coefs[N];
+// for (size_t i = N; i > 0; --i) {
+// bcoef *= x;
+// bcoef += coefs[i - 1];
+// }
+// return bcoef;
+// }
+
+// template <size_t M>
+// PUGS_INLINE constexpr friend void
+// divide(const PolynomialP<N>& P1, const PolynomialP<M>& P2, PolynomialP<N>& Q, PolynomialP<N>& R)
+// {
+// const size_t Nr = P1.realDegree();
+// const size_t Mr = P2.realDegree();
+// R.coefficients() = P1.coefficients();
+// Q.coefficients() = zero;
+// for (ssize_t k = Nr - Mr; k >= 0; --k) {
+// Q.coefficients()[k] = R.coefficients()[Mr + k] / P2.coefficients()[Mr];
+// for (ssize_t j = Mr + k; j >= k; --j) {
+// R.coefficients()[j] -= Q.coefficients()[k] * P2.coefficients()[j - k];
+// }
+// }
+// for (size_t j = Mr; j <= Nr; ++j) {
+// R.coefficients()[j] = 0;
+// }
+// }
+
+// PUGS_INLINE
+// constexpr friend PolynomialP<N + 1>
+// primitive(const PolynomialP<N>& P)
+// {
+// TinyVector<N + 2> coefs;
+// for (size_t i = 0; i < N + 1; ++i) {
+// coefs[i + 1] = P.coefficients()[i] / double(i + 1);
+// }
+// coefs[0] = 0;
+// return PolynomialP<N + 1>{coefs};
+// }
+
+// PUGS_INLINE
+// constexpr friend std::ostream&
+// operator<<(std::ostream& os, const PolynomialP<N>& P)
+// {
+// // os << "P(x) = ";
+// bool all_coef_zero = true;
+// if (N == 0) {
+// os << P.coefficients()[0];
+// return os;
+// }
+// if (N != 1) {
+// if (P.coefficients()[N] != 0.) {
+// if (P.coefficients()[N] < 0.) {
+// os << "- ";
+// }
+// if (P.coefficients()[N] != 1 && P.coefficients()[N] != -1) {
+// os << std::abs(P.coefficients()[N]);
+// }
+// os << "x^" << N;
+// all_coef_zero = false;
+// }
+// }
+// for (size_t i = N - 1; i > 1; --i) {
+// if (P.coefficients()[i] != 0.) {
+// if (P.coefficients()[i] > 0.) {
+// os << " + ";
+// } else if (P.coefficients()[i] < 0.) {
+// os << " - ";
+// }
+// if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) {
+// os << std::abs(P.coefficients()[i]);
+// }
+// os << "x^" << i;
+// all_coef_zero = false;
+// }
+// }
+// if (P.coefficients()[1] != 0.) {
+// if (P.coefficients()[1] > 0. && N != 1) {
+// os << " + ";
+// } else if (P.coefficients()[1] < 0.) {
+// os << " - ";
+// }
+// if (P.coefficients()[1] != 1 && P.coefficients()[1] != -1) {
+// os << std::abs(P.coefficients()[1]);
+// }
+// os << "x";
+// all_coef_zero = false;
+// }
+// if (P.coefficients()[0] != 0. || all_coef_zero) {
+// if (P.coefficients()[0] > 0.) {
+// os << " + ";
+// } else if (P.coefficients()[0] < 0.) {
+// os << " - ";
+// }
+// os << std::abs(P.coefficients()[0]);
+// }
+// return os;
+// }
+
+// PUGS_INLINE
+// constexpr friend void
+// lagrangeBasis(const TinyVector<N + 1> zeros, TinyVector<N + 1, PolynomialP<N>>& basis)
+// {
+// PolynomialP<N> lj;
+// for (size_t j = 0; j < N + 1; ++j) {
+// basis[j] = lagrangePolynomialP(zeros, j);
+// }
+// }
+
+// PUGS_INLINE
+// constexpr friend PolynomialP<N>
+// lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, const TinyVector<N + 1, PolynomialP<N>>& basis)
+// {
+// PolynomialP<N> P(zero);
+// // lagrangeBasis({0, 0, 0}, basis);
+// for (size_t j = 0; j < N + 1; ++j) {
+// P += basis[j] * lagrange_coefs[j];
+// }
+// return P;
+// }
+
+// template <size_t N>
+// PUGS_INLINE constexpr PolynomialP<N> lagrangePolynomialP(const TinyVector<N + 1> zeros, const size_t k);
+
+// template <size_t N>
+// PUGS_INLINE constexpr TinyVector<N, PolynomialP<N - 1>>
+// lagrangeBasis(const TinyVector<N>& zeros)
+// {
+// static_assert(N >= 1, "invalid degree");
+// TinyVector<N, PolynomialP<N - 1>> basis;
+// PolynomialP<N - 1> lj;
+// for (size_t j = 0; j < N; ++j) {
+// basis[j] = lagrangePolynomialP<N - 1>(zeros, j);
+// }
+// return basis;
+// }
+
+// template <size_t N>
+// PUGS_INLINE constexpr double
+// integrate(const PolynomialP<N>& P, const double& xinf, const double& xsup)
+// {
+// PolynomialP<N + 1> Q = primitive(P);
+// return (Q(xsup) - Q(xinf));
+// }
+
+// template <size_t N>
+// PUGS_INLINE constexpr double
+// symmetricIntegrate(const PolynomialP<N>& P, const double& delta)
+// {
+// Assert(delta > 0, "A positive delta is needed for symmetricIntegrate");
+// double integral = 0.;
+// for (size_t i = 0; i <= N; ++i) {
+// if (i % 2 == 0)
+// integral += 2. * P.coefficients()[i] * std::pow(delta, i + 1) / (i + 1);
+// }
+// return integral;
+// }
+
+// template <size_t N>
+// PUGS_INLINE constexpr auto
+// derivative(const PolynomialP<N>& P)
+// {
+// if constexpr (N == 0) {
+// return PolynomialP<0>(0);
+// } else {
+// TinyVector<N> coefs;
+// for (size_t i = 0; i < N; ++i) {
+// coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+// }
+// return PolynomialP<N - 1>(coefs);
+// }
+// }
+
+// template <size_t N>
+// PUGS_INLINE constexpr PolynomialP<N>
+// lagrangePolynomialP(const TinyVector<N + 1> zeros, const size_t k)
+// {
+// for (size_t i = 0; i < N; ++i) {
+// Assert(zeros[i] < zeros[i + 1], "Interpolation values must be strictly increasing in Lagrange polynomialPs");
+// }
+// PolynomialP<N> lk;
+// lk.coefficients() = zero;
+// lk.coefficients()[0] = 1;
+// for (size_t i = 0; i < N + 1; ++i) {
+// if (k == i)
+// continue;
+// double factor = 1. / (zeros[k] - zeros[i]);
+// PolynomialP<1> P({-zeros[i] * factor, factor});
+// lk *= P;
+// }
+// return lk;
+// }
+
+#endif // POLYNOMIALP_HPP
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1f1f2e4964ed417376072f2e33a876996809f947
--- /dev/null
+++ b/tests/test_PolynomialP.cpp
@@ -0,0 +1,211 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <analysis/PolynomialP.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class PolynomialP<0, 2>;
+template class PolynomialP<1, 2>;
+template class PolynomialP<2, 2>;
+template class PolynomialP<3, 2>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PolynomialP", "[analysis]")
+{
+ SECTION("construction")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ REQUIRE_NOTHROW(PolynomialP<2, 2>(coef));
+ }
+
+ SECTION("degree")
+ {
+ TinyVector<3> coef(1, 2, 3);
+ PolynomialP<1, 2> P(coef);
+ REQUIRE(P.degree() == 1);
+ REQUIRE(P.dim() == 2);
+ }
+
+ SECTION("equality")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<3> coef2(1, 2, 3);
+ TinyVector<6> coef3(1, 2, 3, 3, 3, 3);
+ PolynomialP<2, 2> P(coef);
+ PolynomialP<2, 2> Q(coef);
+ PolynomialP<2, 2> R(coef3);
+
+ REQUIRE(P == Q);
+ REQUIRE(P != R);
+ }
+
+ SECTION("addition")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<6> coef2(1, 2, 3, -2, -1, -3);
+ TinyVector<6> coef3(2, 4, 6, 2, 4, 3);
+ PolynomialP<2, 2> P(coef);
+ PolynomialP<2, 2> Q(coef2);
+ PolynomialP<2, 2> R(coef3);
+ REQUIRE(R == (P + Q));
+ REQUIRE((P + Q) == R);
+ }
+
+ SECTION("opposed")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<6> coef2(-1, -2, -3, -4, -5, -6);
+ PolynomialP<2, 2> P(coef);
+ REQUIRE(-P == PolynomialP<2, 2>(coef2));
+ }
+
+ // SECTION("difference")
+ // {
+ // Polynomial<2> P(2, 3, 4);
+ // Polynomial<2> Q(3, 4, 5);
+ // Polynomial<2> D(-1, -1, -1);
+ // REQUIRE(D == (P - Q));
+ // Polynomial<3> R(2, 3, 4, 1);
+ // REQUIRE(D == (P - Q));
+ // REQUIRE((P - R) == Polynomial<3>{0, 0, 0, -1});
+ // R -= P;
+ // REQUIRE(R == Polynomial<3>(0, 0, 0, 1));
+ // }
+
+ // SECTION("product_by_scalar")
+ // {
+ // Polynomial<2> P(2, 3, 4);
+ // Polynomial<2> M(6, 9, 12);
+ // REQUIRE(M == (P * 3));
+ // REQUIRE(M == (3 * P));
+ // }
+
+ // SECTION("product")
+ // {
+ // Polynomial<2> P(2, 3, 4);
+ // Polynomial<3> Q(1, 2, -1, 1);
+ // Polynomial<4> R;
+ // Polynomial<5> S;
+ // R = P;
+ // S = P;
+ // S *= Q;
+ // REQUIRE(Polynomial<5>(2, 7, 8, 7, -1, 4) == (P * Q));
+ // REQUIRE(Polynomial<5>(2, 7, 8, 7, -1, 4) == S);
+ // // REQUIRE_THROWS_AS(R *= Q, AssertError);
+ // }
+
+ // SECTION("divide")
+ // {
+ // Polynomial<2> P(1, 0, 1);
+ // Polynomial<1> Q(0, 1);
+ // Polynomial<1> Q1(0, 1);
+
+ // Polynomial<2> R;
+ // Polynomial<2> S;
+ // REQUIRE(P.realDegree() == 2);
+ // REQUIRE(Q.realDegree() == 1);
+ // REQUIRE(Q1.realDegree() == 1);
+
+ // divide(P, Q1, R, S);
+ // REQUIRE(Polynomial<2>{1, 0, 0} == S);
+ // REQUIRE(Polynomial<2>{0, 1, 0} == R);
+ // }
+
+ // SECTION("evaluation")
+ // {
+ // Polynomial<2> P(2, -3, 4);
+ // REQUIRE(P(3) == 29);
+ // }
+
+ // SECTION("primitive")
+ // {
+ // Polynomial<2> P(2, -3, 4);
+ // TinyVector<4> coefs = zero;
+ // Polynomial<3> Q(coefs);
+ // Q = primitive(P);
+ // Polynomial<3> R(0, 2, -3. / 2, 4. / 3);
+ // REQUIRE(Q == R);
+ // }
+
+ // SECTION("integrate")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // double xinf = -1;
+ // double xsup = 1;
+ // double result = integrate(P, xinf, xsup);
+ // REQUIRE(result == 6);
+ // result = symmetricIntegrate(P, 2);
+ // REQUIRE(result == 24);
+ // }
+
+ // SECTION("derivative")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // Polynomial<1> Q = derivative(P);
+ // REQUIRE(Q == Polynomial<1>(-3, 6));
+
+ // Polynomial<0> P2(3);
+
+ // Polynomial<0> R(0);
+ // REQUIRE(derivative(P2) == R);
+ // }
+
+ // SECTION("affectation")
+ // {
+ // Polynomial<2> Q(2, -3, 3);
+ // Polynomial<4> R(2, -3, 3, 0, 0);
+ // Polynomial<4> P(0, 1, 2, 3, 3);
+ // P = Q;
+ // REQUIRE(P == R);
+ // }
+
+ // SECTION("affectation addition")
+ // {
+ // Polynomial<2> Q(2, -3, 3);
+ // Polynomial<4> R(2, -2, 5, 3, 3);
+ // Polynomial<4> P(0, 1, 2, 3, 3);
+ // P += Q;
+ // REQUIRE(P == R);
+ // }
+
+ // SECTION("power")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // Polynomial<4> R(4, -12, 21, -18, 9);
+ // Polynomial<1> Q(0, 2);
+ // Polynomial<2> S = Q.pow<2>(2);
+ // REQUIRE(P.pow<2>(2) == R);
+ // REQUIRE(S == Polynomial<2>(0, 0, 4));
+ // }
+
+ // SECTION("composition")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // Polynomial<1> Q(0, 2);
+ // Polynomial<2> R(2, -1, 3);
+ // Polynomial<2> S(1, 2, 2);
+ // REQUIRE(P.compose(Q) == Polynomial<2>(2, -6, 12));
+ // REQUIRE(P.compose2(Q) == Polynomial<2>(2, -6, 12));
+ // REQUIRE(R(S) == Polynomial<4>(4, 10, 22, 24, 12));
+ // }
+
+ // SECTION("Lagrange polynomial")
+ // {
+ // Polynomial<1> S(0.5, -0.5);
+ // Polynomial<1> Q;
+ // Q = lagrangePolynomial<1>(TinyVector<2>{-1, 1}, 0);
+ // REQUIRE(S == Q);
+ // Polynomial<2> P(0, -0.5, 0.5);
+ // Polynomial<2> R;
+ // R = lagrangePolynomial<2>(TinyVector<3>{-1, 0, 1}, 0);
+ // REQUIRE(R == P);
+ // const std::array<Polynomial<2>, 3> basis = lagrangeBasis(TinyVector<3>{-1, 0, 1});
+ // REQUIRE(lagrangeToCanonical(TinyVector<3>{1, 0, 1}, basis) == Polynomial<2>(TinyVector<3>{0, 0, 1}));
+ // }
+}