Polynomial.hpp

#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;
  static_assert(N >= 0, "Polynomial degree must be non-negative");

 public:
  PUGS_INLINE
  constexpr size_t
  degree() const
  {
    return N;
  }

  PUGS_INLINE
  constexpr size_t
  realDegree() const
  {
    for (size_t j = N; j > 0; j--) {
      if (std::abs(coefficients()[j]) > 1.e-14) {
        return j;
      }
    }
    return 0;
  }

  PUGS_INLINE
  constexpr const Coefficients&
  coefficients() const
  {
    return m_coefficients;
  }

  PUGS_INLINE
  constexpr Coefficients&
  coefficients()
  {
    return m_coefficients;
  }

  PUGS_INLINE constexpr bool
  operator==(const Polynomial& q) const
  {
    return m_coefficients == q.m_coefficients;
  }

  PUGS_INLINE
  constexpr bool
  operator!=(const Polynomial& q) const
  {
    return not(*this == q);
  }

  template <size_t M>
  PUGS_INLINE constexpr Polynomial<std::max(M, N)>
  operator+(const Polynomial<M>& Q) const
  {
    Polynomial<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;
  }

  PUGS_INLINE constexpr Polynomial
  operator-() const
  {
    return Polynomial{-m_coefficients};
  }

  template <size_t M>
  PUGS_INLINE constexpr Polynomial<std::max(M, N)>
  operator-(const Polynomial<M>& Q) const
  {
    Polynomial<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 Polynomial&
  operator=(const Polynomial<M>& Q)
  {
    coefficients() = zero;
    for (size_t i = N + 1; i <= M; ++i) {
      Assert(Q.coefficients()[i] == 0, "degree of polynomial to small in assignation");
    }
    //    static_assert(N >= M, "degree of polynomial 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 Polynomial&
  operator+=(const Polynomial<M>& Q)
  {
    static_assert(N >= M, "Polynomial 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 Polynomial&
  operator-=(const Polynomial<M>& Q)
  {
    static_assert(N >= M, "Polynomial 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 Polynomial<M + N>
  operator*(const Polynomial<M>& Q) const
  {
    Polynomial<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 Polynomial&
  operator*=(const Polynomial<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 polynomial");
    }
    Polynomial 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 Polynomial
  operator*(const double& lambda) const
  {
    return Polynomial(lambda * m_coefficients);
  }

  template <size_t M>
  PUGS_INLINE constexpr Polynomial<M * N>
  compose(const Polynomial<M>& Q) const
  {
    Polynomial<M * N> P;
    P.coefficients() = zero;
    Polynomial<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 Polynomial<M * N>
  power(const Polynomial<M>& Q) const
  {
    return coefficients()[I] * Q.template pow<N>(I);
  }

  template <size_t M, size_t... I>
  PUGS_INLINE constexpr Polynomial<M * N>
  compose2(const Polynomial<M>& Q, std::index_sequence<I...>) const
  {
    Polynomial<M * N> P;
    P.coefficients() = zero;
    P                = (power<M, I>(Q) + ...);
    return P;
  }

  template <size_t M>
  PUGS_INLINE constexpr Polynomial<M * N>
  compose2(const Polynomial<M>& Q) const
  {
    Polynomial<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 Polynomial<M * N>
  operator()(const Polynomial<M>& Q) const
  {
    Polynomial<M * N> P;
    P.coefficients() = zero;
    Polynomial<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 Polynomial<M * N>
  pow(size_t power) const
  {
    Assert(power <= M, "You declared a polynomial of degree too small for return of the pow function");
    Polynomial<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 Polynomial
  operator*(const double& lambda, const Polynomial& P)
  {
    return P * lambda;
  }

  PUGS_INLINE
  constexpr double
  operator()(double x) const
  {
    double p_x = m_coefficients[N];
    for (size_t i = N; i > 0; --i) {
      p_x *= x;
      p_x += m_coefficients[i - 1];
    }
    return p_x;
  }

  template <size_t M>
  PUGS_INLINE constexpr friend void
  divide(const Polynomial<N>& P1, const Polynomial<M>& P2, Polynomial<N>& Q, Polynomial<N>& R)
  {
    const size_t Nr  = P1.realDegree();
    const size_t Mr  = P2.realDegree();
    R.coefficients() = P1.coefficients();
    Q.coefficients() = zero;
    for (size_t k = Nr - Mr + 1; k > 0; --k) {
      Q.coefficients()[k - 1] = R.coefficients()[Mr + k - 1] / P2.coefficients()[Mr];
      for (size_t j = Mr + k - 1; j > (k - 1); --j) {
        R.coefficients()[j] -= Q.coefficients()[k - 1] * P2.coefficients()[j - k];
      }
    }
    for (size_t j = Mr; j < Nr + 1; ++j) {
      R.coefficients()[j] = 0;
    }
  }

  PUGS_INLINE
  constexpr friend Polynomial<N + 1>
  primitive(const Polynomial& 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& 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, Polynomial<N>>& basis)
  {
    Polynomial<N> lj;
    for (size_t j = 0; j < N + 1; ++j) {
      basis[j] = lagrangePolynomial(zeros, j);
    }
  }

  PUGS_INLINE
  constexpr friend Polynomial<N>
  lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, const std::array<Polynomial<N>, N + 1>& basis)
  {
    Polynomial<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;
  }

  PUGS_INLINE constexpr Polynomial& operator=(const Polynomial&) = default;
  PUGS_INLINE constexpr Polynomial& operator=(Polynomial&&) = default;

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

  PUGS_INLINE
  constexpr Polynomial(const Polynomial&) noexcept = default;

  PUGS_INLINE
  constexpr Polynomial(Polynomial&&) noexcept = default;

  template <typename... T>
  explicit PUGS_INLINE constexpr Polynomial(T&&... coefficients) noexcept : m_coefficients(coefficients...)
  {
    // static_assert(sizeof...(T) == N + 1, "invalid number of parameters");
    // static_assert((std::is_convertible_v<T, double> and ...), "arguments must be convertible to double");
  }

  PUGS_INLINE
  constexpr Polynomial() noexcept = default;

  ~Polynomial() = default;
};

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

template <size_t N>
PUGS_INLINE constexpr Polynomial<N> lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k);

template <size_t N>
PUGS_INLINE constexpr std::array<Polynomial<N - 1>, N>
lagrangeBasis(const TinyVector<N>& zeros)
{
  static_assert(N >= 1, "invalid degree");
  std::array<Polynomial<N - 1>, N> basis;
  for (size_t j = 0; j < N; ++j) {
    basis[j] = lagrangePolynomial<N - 1>(zeros, j);
  }
  return basis;
}

template <size_t N>
PUGS_INLINE constexpr double
integrate(const Polynomial<N>& P, const double& xinf, const double& xsup)
{
  Polynomial<N + 1> Q = primitive(P);
  return (Q(xsup) - Q(xinf));
}

template <size_t N>
PUGS_INLINE constexpr double
symmetricIntegrate(const Polynomial<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 Polynomial<N>& P)
{
  if constexpr (N == 0) {
    return Polynomial<0>(zero);
  } else {
    TinyVector<N> coefs;
    for (size_t i = 0; i < N; ++i) {
      coefs[i] = double(i + 1) * P.coefficients()[i + 1];
    }
    return Polynomial<N - 1>(coefs);
  }
}

template <size_t N>
PUGS_INLINE constexpr Polynomial<N>
lagrangePolynomial(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 polynomials");
  }
  Polynomial<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]);
    Polynomial<1> P(TinyVector<2>{-zeros[i] * factor, factor});
    lk *= P;
  }
  return lk;
}

#endif   // POLYNOMIAL_HPP