Polynomial1D.hpp

#ifndef POLYNOMIAL_1D_HPP
#define POLYNOMIAL_1D_HPP

#include <utils/Array.hpp>
#include <utils/Exceptions.hpp>

#include <algorithm>

class [[nodiscard]] Polynomial1D
{
 private:
  std::vector<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)
  {
    // return std::move(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());
    std::fill(begin(product.m_coefficients), end(product.m_coefficients), 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());
      for (size_t i = 0; i < greater.m_coefficients.size(); ++i) {
        result.m_coefficients[i] = greater.m_coefficients[i];
      }

      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());
    std::fill(begin(ratio.m_coefficients), end(ratio.m_coefficients), 0);
    Polynomial1D remaining(this->degree());
    for (size_t i = 0; i < m_coefficients.size(); ++i) {
      remaining.m_coefficients[i] = m_coefficients[i];
    }

    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];
      }
    }

    if (q_degree == 0) {
      remaining.m_coefficients.resize(1);
      remaining.m_coefficients[0] = 0;
    } else {
      remaining.m_coefficients.resize(std::max(q_degree - 1, 1ul));
    }
    // 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());
    std::fill(begin(ratio.m_coefficients), end(ratio.m_coefficients), 0);
    Polynomial1D remaining(this->degree());
    for (size_t i = 0; i < m_coefficients.size(); ++i) {
      remaining.m_coefficients[i] = m_coefficients[i];
    }

    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(), "invalid coefficient number");
    return m_coefficients[i];
  }

  PUGS_INLINE
  const double& coefficient(const size_t i) const
  {
    Assert(i < m_coefficients.size(), "invalid coefficient number");
    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];
    }

    size_t real_degree = this->_getRealDegree();
    if (real_degree + 1 < m_coefficients.size()) {
      std::vector<double> real_coefficients(real_degree + 1);
      for (size_t i = 0; i <= real_degree; ++i) {
        real_coefficients[i] = m_coefficients[i];
      }
      m_coefficients = std::move(real_coefficients);
    }
  }

  PUGS_INLINE
  explicit Polynomial1D(const size_t degree) : m_coefficients(degree + 1)
  {
    std::fill(m_coefficients.begin(), m_coefficients.end(), 0);
  }

  Polynomial1D(const Polynomial1D&) = default;
  Polynomial1D(Polynomial1D &&)     = default;

  ~Polynomial1D() = default;
};

#endif   // POLYNOMIAL_1D_HPP