Poly.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2026, Davide Stocco and Enrico Bertolazzi.                  *
 *                                                                           *
 * The Sturm project is distributed under the BSD 2-Clause License.          *
 *                                                                           *
 * Davide Stocco                                           Enrico Bertolazzi *
 * University of Trento                                 University of Trento *
 * davide.stocco@unitn.it                         enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef INCLUDE_STURM_POLY_HH
#define INCLUDE_STURM_POLY_HH
 
#include "Sturm.hh"
 
namespace Sturm {

  template <typename Real>
  class Poly : public Eigen::Vector<Real, Eigen::Dynamic> {
   public:
    using Vector =
        Eigen::Vector<Real, Eigen::Dynamic>; 
    constexpr static const Real EPSILON{
      std::numeric_limits<Real>::epsilon()}; 
 
   private:
    Integer m_order;                         

    Vector &to_eigen() {
      return *static_cast<Vector *>(this);
    }
 
   public:
    Poly() : m_order(0) {}

    explicit Poly(Integer order) : m_order(order) {
      this->resize(order);
      this->setZero();
    }

    Poly(const Poly &c) : Vector(c), m_order(c.m_order) {}

    explicit Poly(const Vector &c) {
      this->resize(c.size());
      this->Vector::operator=(c);  // to avoid recursion
      this->m_order = Integer(c.size());
    }

    const Vector &to_eigen() const {
      return *static_cast<const Vector *>(this);
    }

    void set_order(Integer order) {
      this->m_order = order;
      this->resize(order);
      this->setZero();
    }

    void set_degree(Integer degree) {
      this->m_order = degree + 1;
      this->resize(this->m_order);
      this->setZero();
    }

    Poly &set_scalar(Real s) {
      this->resize(1);
      this->coeffRef(0) = s;
      this->m_order     = 1;
      return *this;
    }

    Poly &set_monomial(Real a) {
      this->resize(2);
      this->coeffRef(0) = a;
      this->coeffRef(1) = 1;
      this->m_order     = 2;
      return *this;
    }

    Vector coeffs() const {
      return this->to_eigen();
    }

    Integer degree() const {
      return this->m_order - 1;
    }

    Integer order() const {
      return this->m_order;
    }

    std::string to_string() const {
      if (this->order() <= 0) {
        return "(empty polynomial)";
      }
      if (this->order() == 1) {
        return std::to_string(this->coeff(0));
      };
      if ((*this).cwiseAbs().maxCoeff() == 0) {
        return "0";
      };
 
      bool empty{true};   // true means that all the coefficients are zero
      std::string s{""};  // sign
      Real c{0};          // coefficient
      std::string e{""};  // exponent
      std::string res{""};
 
      // Check if the first coefficient is different from zero
      if (this->coeff(0) != 0) {
        res   = std::to_string(this->coeff(0));
        empty = false;
      }
 
      for (Integer i{1}; i < this->order(); ++i) {
        // If the coefficient is negative...
        if (this->coeff(i) < 0) {
          // and if the previous coefficients are null...
          if (empty) {
            s     = "";  // ...do not write the sign
            c     = this->coeff(i);
            empty = false;
          } else {
            s = " - ";            // ...otherwise write the sign as a separator
            c = -this->coeff(i);  // ...and change the sign of the coefficient
          }
 
        } else if (this->coeff(i) > 0) {
          // If the coefficient is positive...
          c = this->coeff(i);
          // and if the previous coefficients are null...
          if (empty) {
            s     = "";  // ...do not write the sign
            empty = false;
          } else {
            s = " + ";   // ...otherwise write the sign as a separator
          }
 
        } else {
          // If the coefficient is zero, skip to the next
          continue;
        }
 
        // If the exponent is 1, do not write it
        if (i == 1) {
          e = "x";
        } else {
          e = "x^" + std::to_string(i);
        }
 
        // If the coefficient is 1, do not write it
        if (c == 1) {
          res += s;
          res += e;
        } else {
          res += s + std::to_string(c) + e;
        }
      }
      return res;
    }

    Real evaluate(Real x) const {
      Integer n{this->m_order - 1};
      Real res{this->coeff(n)};
      while (n-- > 0) {
        res = res * x + this->coeff(n);
      }
      return res;
    }

    Real evaluate_derivative(Real x) const {
      Integer n{this->m_order - 1};
      Real Dp{this->coeff(n) * n};
      while (--n > 0) {
        Dp = Dp * x + this->coeff(n) * n;
      }
      return Dp;
    }

    void evaluate(Real x, Real &p, Real &Dp) const {
      Integer n{this->m_order - 1};
      p  = this->coeff(n);
      Dp = this->coeff(n) * n;
      while (--n > 0) {
        p  = p * x + this->coeff(n);
        Dp = Dp * x + this->coeff(n) * n;
      }
      p = p * x + this->coeff(0);
    }

    Real leading_coeff() const {
      return this->coeff(m_order - 1);
    }

    void derivative(Poly &res) const {
      res.resize(this->m_order - 1);
      for (Integer i{1}; i < this->m_order; ++i) {
        res.coeffRef(i - 1) = i * this->coeff(i);
      }
      res.m_order = this->m_order - 1;
    }

    void integral(Poly &res) const {
      res.resize(this->m_order + 1);
      res.coeffRef(0) = 0;
      for (Integer i{1}; i <= this->m_order; ++i) {
        res.coeffRef(i) = this->coeff(i - 1) / i;
      }
      res.m_order = this->m_order + 1;
    }

    void integral(Poly &res, Real c) const {
      res.resize(this->m_order + 1);
      res.coeffRef(0) = c;
      for (Integer i{1}; i <= this->m_order; ++i) {
        res.coeffRef(i) = this->coeff(i - 1) / i;
      }
      res.m_order = this->m_order + 1;
    }

    Real normalize() {
      Real scale{this->cwiseAbs().maxCoeff()};
      if (scale > static_cast<Real>(0.0)) {
        this->to_eigen() /= scale;
      }
      return scale;
    }

    void purge(const Real eps) {
      if (this->m_order > 0) {
        Real max{this->cwiseAbs().maxCoeff()};
        if (max < 1) {
          max = 1;
        }
        Real eps_tmp{eps * max};
        for (Integer i{0}; i < this->m_order; ++i) {
          Real &c_i{this->coeffRef(i)};
          if (std::abs(c_i) <= eps_tmp) {
            c_i = 0;
          }
        }
      }
      this->adjust_degree();
    }

    void adjust_degree() {
      while (this->m_order > 0 && this->coeff(this->m_order - 1) == 0) {
        --this->m_order;
      }
      this->conservativeResize(this->m_order);
    }

    Integer sign_variations() const {
      Integer sign_var{0};
      Integer last_sign{0};
      for (Integer i{0}; i < this->m_order; ++i) {
        Real v = this->coeff(i);
        if (v > 0) {
          if (last_sign == -1) {
            ++sign_var;
          }
          last_sign = 1;
        } else if (v < 0) {
          if (last_sign == 1) {
            ++sign_var;
          }
          last_sign = -1;
        }
      }
      return sign_var;
    }

    void make_monic() {
      this->to_eigen() /= this->coeff(this->m_order - 1);
      this->coeffRef(this->m_order - 1) = 1;
    }

    Poly &operator=(const Poly &p) {
      this->resize(p.m_order);
      this->to_eigen().noalias() = p.to_eigen();
      this->m_order              = p.m_order;
      return *this;
    }

    Poly operator-() {
      return Poly(-this->to_eigen());
    }

    Poly &operator+=(const Poly &p) {
      Integer max_order{std::max(this->m_order, p.m_order)};
      Integer min_order{std::min(this->m_order, p.m_order)};
      this->conservativeResize(
          max_order);  // Resize the vector without destroying the content
      this->head(min_order).noalias() += p.head(min_order);
      Integer n_tail{p.m_order - this->m_order};
      if (n_tail > 0) {
        this->tail(n_tail).noalias() = p.tail(n_tail);
      }
      this->m_order = max_order;
      return *this;
    }

    Poly &operator-=(const Poly &p) {
      Integer max_order{std::max(this->m_order, p.m_order)};
      Integer min_order{std::min(this->m_order, p.m_order)};
      this->conservativeResize(
          max_order);  // Resize the vector without destroying the content
      this->head(min_order).noalias() -= p.head(min_order);
      Integer n_tail{p.m_order - this->m_order};
      if (n_tail > 0) {
        this->tail(n_tail).noalias() = -p.tail(n_tail);
      }
      this->m_order = max_order;
      return *this;
    }

    Poly &operator*=(const Poly &p) {
      Vector a(this->to_eigen());
      Integer new_order{this->m_order + p.m_order - 1};
      this->resize(new_order);
      this->setZero();
      for (Integer i{0}; i < this->m_order; ++i) {
        for (Integer j{0}; j < p.m_order; ++j) {
          this->coeffRef(i + j) += a.coeff(i) * p.coeff(j);
        }
      }
      this->m_order = new_order;
      return *this;
    }

    Poly &operator+=(Real s) {
      if (this->m_order > 0)
        this->coeffRef(0) += s;
      else {
        this->resize(1);
        this->coeffRef(0) = s;
        this->m_order     = 1;
      }
      return *this;
    }

    Poly &operator-=(Real s) {
      if (this->m_order > 0) {
        this->coeffRef(0) -= s;
      } else {
        this->resize(1);
        this->coeffRef(0) = -s;
        this->m_order     = 1;
      }
      return *this;
    }

    Poly &operator*=(Real s) {
      this->to_eigen() *= s;
      return *this;
    }
 
  };  // class Poly

  template <typename Real>
  inline Poly<Real> operator+(const Poly<Real> &p_1, const Poly<Real> &p_2) {
    Integer max_order{std::max(p_1.order(), p_2.order())};
    Integer min_order{std::min(p_1.order(), p_2.order())};
    Poly<Real> res(max_order);
    res.head(min_order).noalias() = p_1.head(min_order) + p_2.head(min_order);
    Integer n_tail{max_order - min_order};
    if (n_tail > 0) {
      if (p_1.order() > p_2.order()) {
        res.tail(n_tail).noalias() = p_1.tail(n_tail);
      } else {
        res.tail(n_tail).noalias() = p_2.tail(n_tail);
      }
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator+(const Poly<Real> &p, Real s) {
    Integer max_order{std::max(p.order(), 1)};
    Poly<Real> res(max_order);
    if (p.order() > 0) {
      res.coeffRef(0) = p.coeff(0) + s;
      if (p.order() > 1)
        res.tail(p.order() - 1).noalias() = p.tail(p.order() - 1);
    } else {
      res.coeffRef(0) = s;
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator+(Real s, const Poly<Real> &p) {
    Integer max_order{std::max(p.order(), 1)};
    Poly<Real> res(max_order);
    if (p.order() > 0) {
      res.coeffRef(0) = s + p.coeff(0);
      if (p.order() > 1)
        res.tail(p.order() - 1).noalias() = p.tail(p.order() - 1);
    } else {
      res.coeffRef(0) = s;
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator-(const Poly<Real> &p_1, const Poly<Real> &p_2) {
    Integer max_order{std::max(p_1.order(), p_2.order())};
    Integer min_order{std::min(p_1.order(), p_2.order())};
    Poly<Real> res(max_order);
    res.head(min_order).noalias() = p_1.head(min_order) - p_2.head(min_order);
    Integer n_tail{max_order - min_order};
    if (n_tail > 0) {
      if (p_1.order() > p_2.order()) {
        res.tail(n_tail).noalias() = p_1.tail(n_tail);
      } else {
        res.tail(n_tail).noalias() = -p_2.tail(n_tail);
      }
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator-(const Poly<Real> &p, Real s) {
    Integer max_order{std::max(p.order(), 1)};
    Poly<Real> res(max_order);
    if (p.order() > 0) {
      res.coeffRef(0) = p.coeff(0) - s;
      if (p.order() > 1)
        res.tail(p.order() - 1).noalias() = p.tail(p.order() - 1);
    } else {
      res.coeffRef(0) = -s;
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator-(Real s, const Poly<Real> &p) {
    Integer max_order{std::max(p.order(), 1)};
    Poly<Real> res(max_order);
    if (p.order() > 0) {
      res.coeffRef(0) = s - p.coeff(0);
      if (p.order() > 1) {
        res.tail(p.order() - 1).noalias() = -p.tail(p.order() - 1);
      }
    } else {
      res.coeffRef(0) = s;
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator*(const Poly<Real> &p_1, const Poly<Real> &p_2) {
    Poly<Real> res(p_1.order() + p_2.order() - 1);
    for (Integer i{0}; i < p_1.order(); ++i) {
      for (Integer j{0}; j < p_2.order(); ++j) {
        res.coeffRef(i + j) += p_1.coeff(i) * p_2.coeff(j);
      }
    }
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator*(Real p, const Poly<Real> &s) {
    Poly<Real> res(s.order());
    res.noalias() = p * s.to_eigen();
    return res;
  }

  template <typename Real>
  inline Poly<Real> operator*(const Poly<Real> &s, Real p) {
    Poly<Real> res(s.order());
    res.noalias() = s.to_eigen() * p;
    return res;
  }

  template <typename Real>
  inline void divide(const Poly<Real> &p_1,
                     const Poly<Real> &p_2,
                     Poly<Real> &q,
                     Poly<Real> &r) {
    // Scale polynomials as p_1_norm(x) = p_1(x) / scale_p_1, and p_2_norm(x) =
    // p_2(x) / scale_p_2
    Poly<Real> p_1_norm(p_1), p_2_norm(p_2);
    Real scale_p_1{p_1_norm.normalize()};
    Real scale_p_2{p_2_norm.normalize()};
 
    // p_1_norm(x) = p_2_norm(x) * q(x) + r(x)
    r = p_1_norm;
    Real leading_b_norm{p_2_norm.leading_coeff()};
    Integer d{r.order() - p_2_norm.order()};
    if (d < 0) {
      // p_1_norm(x) = p_2_norm(x) + r(x)
      q.set_scalar(1);
      r = p_1_norm - p_2_norm;
    } else {
      Integer r_degree{r.degree()};
      q.set_order(d + 1);
 
      STURM_ASSERT(
          std::abs(leading_b_norm) > Poly<Real>::EPSILON,
          "Sturm::Poly::divide(...): leading coefficient of p_2(x) is 0.");
 
      while (d >= 0 && r_degree >= 0) {
        Real leading_r = {r(r_degree)};
        Real ratio{leading_r / leading_b_norm};
        q.coeffRef(d) = ratio;
        r.segment(d, p_2_norm.degree()).noalias() -=
            ratio * p_2_norm.head(p_2_norm.degree());
        r.coeffRef(r_degree) = 0;
        --r_degree;
        --d;
      }
 
      // Do not purge remainder: this can be done externally with r.purge(eps);
      r.adjust_degree();
    }
 
    // Scale back polinomials:
    // - p_1_norm(x) = p_2_norm(x) * q(x) + r(x)
    // - p_1(x) / scale_p_1 = p_2(x) / scale_p_2 * q(x) + r(x)
    // - p_1(x) = p_2(x) * (scale_p_1/scale_p_2) * q(x) + scale_p_1*r(x)
    q *= scale_p_1 / scale_p_2;
    r *= scale_p_1;
  }

  template <typename Real>
  inline void GCD(const Poly<Real> &p_1,
                  const Poly<Real> &p_2,
                  Poly<Real> &gcd,
                  Real eps = Poly<Real>::EPSILON) {
    if (p_2.order() > 0) {
      Poly<Real> q, r;
      Sturm::divide(p_1, p_2, q, r);
      r.purge(eps);
      Sturm::GCD(p_2, r, gcd, eps);
    } else {
      gcd = p_1;
    }
    gcd.normalize();
  }

  template <typename Real>
  inline std::ostream &operator<<(std::ostream &os, const Poly<Real> &p) {
    os << p.to_string();
    return os;
  }
 
}  // namespace Sturm
 
#endif  // INCLUDE_STURM_POLY_HH