diff --git a/src/analysis/Polynomial1D.hpp b/src/analysis/Polynomial1D.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..20d4cf90ac11c8b0ae6c264b9f3ccf632d0c493e
--- /dev/null
+++ b/src/analysis/Polynomial1D.hpp
@@ -0,0 +1,280 @@
+#ifndef POLYNOMIAL_1D_HPP
+#define POLYNOMIAL_1D_HPP
+
+#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
+
+class [[nodiscard]] Polynomial1D
+{
+ private:
+ Array<double> m_coefficients;
+
+ PUGS_INLINE
+ size_t _getRealDegree() const
+ {
+ size_t real_degree = this->degree();
+ while (real_degree > 0 and std::abs(m_coefficients[real_degree]) < 1E-14) {
+ --real_degree;
+ }
+ return real_degree;
+ }
+
+ PUGS_INLINE
+ friend Polynomial1D _simplify(Polynomial1D && p)
+ {
+ size_t real_degree = p._getRealDegree();
+
+ if (real_degree != p.degree()) {
+ Polynomial1D q(real_degree);
+ for (size_t i = 0; i <= real_degree; ++i) {
+ q.coefficient(i) = p.coefficient(i);
+ }
+ return q;
+ } else {
+ return std::move(p);
+ }
+ }
+
+ public:
+ friend std::ostream& operator<<(std::ostream& os, const Polynomial1D& p)
+ {
+ bool has_written = false;
+ for (size_t i = 0; i < p.m_coefficients.size(); ++i) {
+ const double& coef = p.m_coefficients[i];
+ if (coef != 0) {
+ if (coef > 0) {
+ if (has_written) {
+ os << " + ";
+ }
+ os << coef;
+ } else {
+ if (has_written) {
+ os << " - " << -coef;
+ } else {
+ os << coef;
+ }
+ }
+ if (i > 0) {
+ os << "*x";
+ if (i > 1) {
+ os << '^' << i;
+ }
+ }
+ has_written = true;
+ }
+ }
+ if (not has_written) {
+ os << 0;
+ }
+ return os;
+ }
+
+ PUGS_INLINE
+ Polynomial1D operator-() const
+ {
+ Polynomial1D opposite(this->degree());
+ for (size_t i = 0; i < m_coefficients.size(); ++i) {
+ opposite.m_coefficients[i] = -m_coefficients[i];
+ }
+ return _simplify(std::move(opposite));
+ }
+
+ PUGS_INLINE
+ friend Polynomial1D operator*(const double a, const Polynomial1D& p)
+ {
+ Polynomial1D product(p.degree());
+ for (size_t i = 0; i < p.m_coefficients.size(); ++i) {
+ product.m_coefficients[i] = a * p.m_coefficients[i];
+ }
+ return _simplify(std::move(product));
+ }
+
+ PUGS_INLINE
+ Polynomial1D operator*(const Polynomial1D& p) const
+ {
+ Polynomial1D product(this->degree() + p.degree());
+ product.m_coefficients.fill(0);
+ for (size_t i = 0; i < this->m_coefficients.size(); ++i) {
+ for (size_t j = 0; j < p.m_coefficients.size(); ++j) {
+ product.m_coefficients[i + j] += this->m_coefficients[i] * p.m_coefficients[j];
+ }
+ }
+ return _simplify(std::move(product));
+ }
+
+ PUGS_INLINE
+ Polynomial1D operator+(const Polynomial1D& p) const
+ {
+ auto sum_left_is_bigger = [](const Polynomial1D& greater, const Polynomial1D& smaller) {
+ Polynomial1D result(greater.degree());
+ copy_to(greater.m_coefficients, result.m_coefficients);
+ for (size_t i = 0; i < smaller.m_coefficients.size(); ++i) {
+ result.m_coefficients[i] += smaller.m_coefficients[i];
+ }
+ return _simplify(std::move(result));
+ };
+
+ if (m_coefficients.size() >= p.m_coefficients.size()) {
+ return sum_left_is_bigger(*this, p);
+ } else {
+ return sum_left_is_bigger(p, *this);
+ }
+ }
+
+ PUGS_INLINE
+ Polynomial1D operator-(const Polynomial1D& p) const
+ {
+ if (m_coefficients.size() >= p.m_coefficients.size()) {
+ Polynomial1D result(this->degree());
+ for (size_t i = 0; i < p.m_coefficients.size(); ++i) {
+ result.m_coefficients[i] = m_coefficients[i] - p.m_coefficients[i];
+ }
+ for (size_t i = p.m_coefficients.size(); i < m_coefficients.size(); ++i) {
+ result.m_coefficients[i] = m_coefficients[i];
+ }
+ return _simplify(std::move(result));
+ } else {
+ Polynomial1D result(p.degree());
+ for (size_t i = 0; i < m_coefficients.size(); ++i) {
+ result.m_coefficients[i] = m_coefficients[i] - p.m_coefficients[i];
+ }
+ for (size_t i = m_coefficients.size(); i < p.m_coefficients.size(); ++i) {
+ result.m_coefficients[i] = -p.m_coefficients[i];
+ }
+ return _simplify(std::move(result));
+ }
+ }
+
+ PUGS_INLINE
+ Polynomial1D operator%(const Polynomial1D& q) const
+ {
+ const Polynomial1D& p = *this;
+ Polynomial1D ratio(this->degree());
+ ratio.m_coefficients.fill(0);
+ Polynomial1D remaining(this->degree());
+ copy_to(m_coefficients, remaining.m_coefficients);
+
+ const size_t p_degree = p.degree();
+ const size_t q_degree = q._getRealDegree();
+
+ for (ssize_t i = p_degree - q_degree; i >= 0; --i) {
+ ratio.m_coefficients[i] = remaining.m_coefficients[q_degree + i] / q.m_coefficients[q_degree];
+ for (ssize_t j = q_degree + i; j >= i; --j) {
+ remaining.m_coefficients[j] -= ratio.m_coefficients[i] * q.m_coefficients[j - i];
+ }
+ }
+
+ for (size_t j = q_degree; j <= remaining.degree(); ++j) {
+ remaining.m_coefficients[j] = 0;
+ }
+
+ return _simplify(std::move(remaining));
+ }
+
+ PUGS_INLINE
+ Polynomial1D operator/(const Polynomial1D& q) const
+ {
+ const Polynomial1D& p = *this;
+ Polynomial1D ratio(this->degree());
+ ratio.m_coefficients.fill(0);
+ Polynomial1D remaining(this->degree());
+ copy_to(m_coefficients, remaining.m_coefficients);
+
+ const size_t p_degree = p.degree();
+ const size_t q_degree = q._getRealDegree();
+
+ for (ssize_t i = p_degree - q_degree; i >= 0; --i) {
+ ratio.m_coefficients[i] = remaining.m_coefficients[q_degree + i] / q.m_coefficients[q_degree];
+ for (ssize_t j = q_degree + i; j >= i; --j) {
+ remaining.m_coefficients[j] -= ratio.m_coefficients[i] * q.m_coefficients[j - i];
+ }
+ }
+
+ return _simplify(std::move(ratio));
+ }
+
+ PUGS_INLINE size_t degree() const
+ {
+ Assert(m_coefficients.size() > 0);
+ return m_coefficients.size() - 1;
+ }
+
+ PUGS_INLINE
+ double& coefficient(const size_t i)
+ {
+ Assert(i < m_coefficients.size());
+ return m_coefficients[i];
+ }
+
+ PUGS_INLINE
+ const double& coefficient(const size_t i) const
+ {
+ Assert(i < m_coefficients.size());
+ return m_coefficients[i];
+ }
+
+ PUGS_INLINE
+ friend Polynomial1D derive(const Polynomial1D& p)
+ {
+ if (p.degree() > 0) {
+ Polynomial1D derivative{p.degree() - 1};
+ for (size_t i = 0; i <= derivative.degree(); ++i) {
+ derivative.m_coefficients[i] = (i + 1) * p.m_coefficients[i + 1];
+ }
+ return derivative;
+ } else {
+ return Polynomial1D(std::vector<double>{0});
+ }
+ }
+
+ PUGS_INLINE
+ friend Polynomial1D primitive(const Polynomial1D& p)
+ {
+ if (p.degree() > 0) {
+ Polynomial1D primitive{p.degree() + 1};
+ primitive.m_coefficients[0] = 0;
+ for (size_t i = 0; i <= p.degree(); ++i) {
+ primitive.m_coefficients[i + 1] = p.m_coefficients[i] / (i + 1);
+ }
+ return _simplify(std::move(primitive));
+ } else {
+ return Polynomial1D(std::vector<double>{0});
+ }
+ }
+
+ PUGS_INLINE
+ double operator()(const double x) const
+ {
+ double value = m_coefficients[this->degree()];
+ for (ssize_t i = this->degree() - 1; i >= 0; --i) {
+ value *= x;
+ value += m_coefficients[i];
+ }
+ return value;
+ }
+
+ Polynomial1D& operator=(Polynomial1D&&) = default;
+ Polynomial1D& operator=(const Polynomial1D&) = default;
+
+ PUGS_INLINE
+ explicit Polynomial1D(const std::vector<double>& coeffients) : m_coefficients(coeffients.size())
+ {
+ Assert(coeffients.size() > 0, "empty list is not allowed");
+ for (size_t i = 0; i < coeffients.size(); ++i) {
+ m_coefficients[i] = coeffients[i];
+ }
+ }
+
+ PUGS_INLINE
+ explicit Polynomial1D(const size_t degree) : m_coefficients(degree + 1)
+ {
+ m_coefficients.fill(0);
+ }
+
+ Polynomial1D(const Polynomial1D&) = default;
+ Polynomial1D(Polynomial1D &&) = default;
+
+ ~Polynomial1D() = default;
+};
+
+#endif // POLYNOMIAL_1D_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4bb227aac16a433471caa45bdfb8b8bc48f4f2bb..827da322463b0b94ef9ab0ef1aa6714211c53c27 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -83,6 +83,7 @@ add_executable (unit_tests
test_NameProcessor.cpp
test_OStreamProcessor.cpp
test_ParseError.cpp
+ test_Polynomial1D.cpp
test_PugsAssert.cpp
test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp
diff --git a/tests/test_Polynomial1D.cpp b/tests/test_Polynomial1D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..110cfadf74bad94ce8a1e8fc3ff6dff8247cad3f
--- /dev/null
+++ b/tests/test_Polynomial1D.cpp
@@ -0,0 +1,259 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/Polynomial1D.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Polynomial1D", "[analysis]")
+{
+ SECTION("Constructor")
+ {
+ Polynomial1D p({-1, 0, 0, 1});
+
+ REQUIRE(p.coefficient(0) == -1);
+ REQUIRE(p.coefficient(1) == 0);
+ REQUIRE(p.coefficient(2) == 0);
+ REQUIRE(p.coefficient(3) == 1);
+
+ REQUIRE(p.degree() == 3);
+ }
+
+ SECTION("value")
+ {
+ Polynomial1D p({-1, -1.3, 2.7, 1});
+
+ REQUIRE(p.degree() == 3);
+
+ double x = 1.3;
+ REQUIRE(p(x) == Catch::Approx(-1 - 1.3 * x + 2.7 * x * x + x * x * x));
+
+ x = 0;
+ REQUIRE(p(x) == Catch::Approx(-1 - 1.3 * x + 2.7 * x * x + x * x * x));
+
+ x = 7;
+ REQUIRE(p(x) == Catch::Approx(-1 - 1.3 * x + 2.7 * x * x + x * x * x));
+
+ x = -3;
+ REQUIRE(p(x) == Catch::Approx(-1 - 1.3 * x + 2.7 * x * x + x * x * x));
+ }
+
+ SECTION("opposite")
+ {
+ Polynomial1D p({-1, -1.3, 2.7, 1});
+ Polynomial1D r = -p;
+
+ REQUIRE(r.degree() == 3);
+
+ REQUIRE(r.coefficient(0) == 1);
+ REQUIRE(r.coefficient(1) == 1.3);
+ REQUIRE(r.coefficient(2) == -2.7);
+ REQUIRE(r.coefficient(3) == -1);
+ }
+
+ SECTION("sum")
+ {
+ Polynomial1D p({-1, 0, 0, 1});
+ Polynomial1D q({-1, 2, 1.4, 2, -1});
+
+ Polynomial1D r = p + q;
+ REQUIRE(r.degree() == 4);
+
+ REQUIRE(r.coefficient(0) == -2);
+ REQUIRE(r.coefficient(1) == 2);
+ REQUIRE(r.coefficient(2) == 1.4);
+ REQUIRE(r.coefficient(3) == 3);
+ REQUIRE(r.coefficient(4) == -1);
+ }
+
+ SECTION("difference")
+ {
+ Polynomial1D p({-1, 0, 0, 1});
+ Polynomial1D q({-1, 2, 1.4, 2, -1});
+ {
+ Polynomial1D r = p - q;
+ REQUIRE(r.degree() == 4);
+
+ REQUIRE(r.coefficient(0) == 0);
+ REQUIRE(r.coefficient(1) == -2);
+ REQUIRE(r.coefficient(2) == -1.4);
+ REQUIRE(r.coefficient(3) == -1);
+ REQUIRE(r.coefficient(4) == 1);
+ }
+
+ {
+ Polynomial1D s = q - p;
+ REQUIRE(s.degree() == 4);
+
+ REQUIRE(s.coefficient(0) == 0);
+ REQUIRE(s.coefficient(1) == 2);
+ REQUIRE(s.coefficient(2) == 1.4);
+ REQUIRE(s.coefficient(3) == 1);
+ REQUIRE(s.coefficient(4) == -1);
+ }
+
+ {
+ Polynomial1D s = q - q;
+ REQUIRE(s.degree() == 0);
+
+ REQUIRE(s.coefficient(0) == 0);
+ }
+ }
+
+ SECTION("product")
+ {
+ Polynomial1D p({-1, 1, 0, 2});
+ Polynomial1D q({-1, 2, 1.4, 2, -1});
+
+ {
+ Polynomial1D s = 2 * p;
+ REQUIRE(s.degree() == 3);
+
+ REQUIRE(s.coefficient(0) == -2);
+ REQUIRE(s.coefficient(1) == 2);
+ REQUIRE(s.coefficient(2) == 0);
+ REQUIRE(s.coefficient(3) == 4);
+ }
+
+ {
+ Polynomial1D r = p * q;
+ REQUIRE(r.degree() == 7);
+
+ REQUIRE(r.coefficient(0) == Catch::Approx(1));
+ REQUIRE(r.coefficient(1) == Catch::Approx(-3));
+ REQUIRE(r.coefficient(2) == Catch::Approx(0.6));
+ REQUIRE(r.coefficient(3) == Catch::Approx(-2.6));
+ REQUIRE(r.coefficient(4) == Catch::Approx(7));
+ REQUIRE(r.coefficient(5) == Catch::Approx(1.8));
+ REQUIRE(r.coefficient(6) == Catch::Approx(4));
+ REQUIRE(r.coefficient(7) == Catch::Approx(-2));
+ }
+ }
+
+ SECTION("divide")
+ {
+ SECTION("exact")
+ {
+ Polynomial1D p({-1, 1, 0, 2});
+ Polynomial1D q({1, -3, 0.6, -2.6, 7, 1.8, 4, -2});
+
+ {
+ Polynomial1D s = q / p;
+
+ REQUIRE(s.degree() == 4);
+
+ REQUIRE(s.coefficient(0) == -1);
+ REQUIRE(s.coefficient(1) == 2);
+ REQUIRE(s.coefficient(2) == 1.4);
+ REQUIRE(s.coefficient(3) == 2);
+ REQUIRE(s.coefficient(4) == -1);
+ }
+ }
+
+ SECTION("with remainder")
+ {
+ Polynomial1D p({-1, 0, 1});
+ Polynomial1D q({0, 1});
+
+ {
+ Polynomial1D s = p / q;
+
+ REQUIRE(s.degree() == 1);
+
+ REQUIRE(s.coefficient(0) == 0);
+ REQUIRE(s.coefficient(1) == 1);
+ }
+ }
+ }
+
+ SECTION("mod")
+ {
+ SECTION("exact")
+ {
+ Polynomial1D p({-1, 1, 0, 2});
+ Polynomial1D q({1, -3, 0.6, -2.6, 7, 1.8, 4, -2});
+
+ {
+ Polynomial1D s = q % p;
+
+ REQUIRE(s.degree() == 0);
+ REQUIRE(s.coefficient(0) == 0);
+ }
+ }
+
+ SECTION("with remainder")
+ {
+ Polynomial1D p({-1, 0, 1});
+ Polynomial1D q({0, 1});
+
+ {
+ Polynomial1D s = p % q;
+
+ REQUIRE(s.degree() == 0);
+ REQUIRE(s.coefficient(0) == -1);
+ }
+ }
+ }
+
+ SECTION("derive")
+ {
+ SECTION("constant")
+ {
+ Polynomial1D p(std::vector<double>{-1});
+ Polynomial1D p_prime = derive(p);
+
+ REQUIRE(p_prime.degree() == 0);
+ REQUIRE(p_prime.coefficient(0) == 0);
+ }
+
+ SECTION("non constant")
+ {
+ Polynomial1D p({-1, 1, 0, 2});
+ Polynomial1D p_prime = derive(p);
+
+ REQUIRE(p_prime.degree() == 2);
+ REQUIRE(p_prime.coefficient(0) == 1);
+ REQUIRE(p_prime.coefficient(1) == 0);
+ REQUIRE(p_prime.coefficient(2) == 6);
+ }
+ }
+
+ SECTION("primitive")
+ {
+ SECTION("zero")
+ {
+ Polynomial1D p(std::vector<double>{0});
+ Polynomial1D P = primitive(p);
+
+ REQUIRE(P.degree() == 0);
+ REQUIRE(P.coefficient(0) == 0);
+ }
+
+ SECTION("non zero")
+ {
+ Polynomial1D p({-1, 1, 0, 2});
+ Polynomial1D P = primitive(p);
+
+ REQUIRE(P.degree() == 4);
+ REQUIRE(P.coefficient(0) == 0);
+ REQUIRE(P.coefficient(1) == -1);
+ REQUIRE(P.coefficient(2) == 0.5);
+ REQUIRE(P.coefficient(3) == 0);
+ REQUIRE(P.coefficient(4) == 0.5);
+ }
+ }
+
+ SECTION("Output")
+ {
+ Polynomial1D p({-1, 0, 0, 1});
+
+ std::ostringstream os;
+ os << p;
+
+ REQUIRE(os.str() == "-1 + 1*x^3");
+ }
+
+#ifndef NDEBUG
+ SECTION("checking for bounds validation") {}
+#endif // NDEBUG
+}