#ifndef POLYNOMIAL_HPP
#define POLYNOMIAL_HPP

#include <algebra/TinyVector.hpp>

template <size_t N>
class Polynomial
{
 private:
  using Coefficients = TinyVector<N + 1, double>;
  Coefficients m_coefficients;

 public:
  PUGS_INLINE
  constexpr const TinyVector<N + 1>&
  coefficients() const
  {
    return m_coefficients;
  }

  PUGS_INLINE
  constexpr TinyVector<N + 1>&
  coefficients()
  {
    return m_coefficients;
  }

  PUGS_INLINE
  constexpr bool
  operator==(const Polynomial<N>& q) const
  {
    if (m_coefficients != q.coefficients()) {
      return false;
    }
    return true;
  }

  PUGS_INLINE
  constexpr bool
  operator!=(const Polynomial<N>& q) const
  {
    return not this->operator==(q);
  }

  PUGS_INLINE
  constexpr Polynomial<N>
  operator+(const Polynomial<N>& q) const
  {
    TinyVector<N + 1> sum_coefs = m_coefficients + q.coefficients();
    Polynomial<N> S(sum_coefs);
    return S;
  }

  PUGS_INLINE
  constexpr Polynomial<N>
  operator-(const Polynomial<N>& q) const
  {
    TinyVector<N + 1> diff_coefs = m_coefficients - q.coefficients();
    Polynomial<N> D(diff_coefs);
    return D;
  }

  PUGS_INLINE
  constexpr Polynomial<N> operator*(const double& lambda) const
  {
    TinyVector<N + 1> mult_coefs = lambda * m_coefficients;
    Polynomial<N> M(mult_coefs);
    return M;
  }

  PUGS_INLINE
  constexpr friend Polynomial<N> operator*(const double& lambda, const Polynomial<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;
  }

  PUGS_INLINE
  constexpr friend Polynomial<N + 1>
  primitive(const Polynomial<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 Polynomial<N + 1>{coefs};
  }

  PUGS_INLINE
  constexpr friend std::ostream&
  operator<<(std::ostream& os, const Polynomial<N>& P)
  {
    //    os << "P(x) = ";
    bool all_coef_zero = true;
    for (size_t i = N; i > 1; --i) {
      if (P.coefficients()[i] != 0.) {
        if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) {
          os << std::abs(P.coefficients()[i]);
        }
        os << "x^" << i;
        if (P.coefficients()[i - 1] > 0.)
          os << " + ";
        else
          os << " - ";
        all_coef_zero = false;
      } else {
        os << " ";
      }
    }
    if (P.coefficients()[1] != 0.) {
      if (P.coefficients()[1] != 1 && P.coefficients()[1] != -1) {
        os << std::abs(P.coefficients()[1]);
      }
      os << "x";
      if (P.coefficients()[0] > 0.)
        os << " + ";
      else if (P.coefficients()[0] < 0.)
        os << " - ";
      all_coef_zero = false;
    } else {
      os << " ";
    }
    if (P.coefficients()[0] != 0. || all_coef_zero)
      os << std::abs(P.coefficients()[0]);
    return os;
  }

  PUGS_INLINE
  constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {}

  PUGS_INLINE
  constexpr Polynomial(TinyVector<N + 1>&& coefficients) noexcept : m_coefficients{coefficients} {}

  PUGS_INLINE
  constexpr Polynomial() noexcept = default;
  ~Polynomial()                   = default;
};

#endif   // POLYNOMIAL_HPP