Add polynomial 1d class

src/analysis/Polynomial1D.hpp

+#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

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

tests/test_Polynomial1D.cpp

+#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
+}