polynomial.hpp

#ifndef __POLYNOMIAL_H_
#define __POLYNOMIAL_H_

#include <initializer_list>
using std::initializer_list;
#include <vector>
using std::vector;

#include <iostream>
using std::ostream;

#include "rings.hpp"

namespace Rings {

typedef int32_t DEGREE_TYPE;

template<typename R>
class Polynomial : virtual public Commutative_Ring_Element {

protected:

  DEGREE_TYPE deg; 
  vector<R> coeffs; 
  static Polynomial<R> last_divisor;
  static Polynomial<R> last_quotient;
  static Polynomial<R> last_remainder;

  void clear_coeffs();
  void set_degree(DEGREE_TYPE);
  void verify_degree();
  void common_division(const Polynomial<R> &) const;

public:

  // constructors

  Polynomial();
  Polynomial(const Polynomial<R> &);
  Polynomial(const initializer_list<R> &);
  Polynomial(const vector<R> &);

  // required methods

  DEGREE_TYPE degree() const;
  const R & coeff(DEGREE_TYPE) const;
  const R & operator () (const R &) const;

  void set_coeff(DEGREE_TYPE, const R &);

  const Polynomial<R> & operator = (const Ring_Element &);

  virtual bool operator == (const Ring_Element &) const override;
  virtual bool operator != (const Ring_Element &) const override;
  bool operator == (const R &) const;

  const Polynomial<R> & operator / (const Polynomial<R> &) const;
  const Polynomial<R> & operator % (const Polynomial<R> &) const;

  // overriding Commutative_Ring_Element pure virtual

  virtual bool is_one() const override;
  virtual bool is_zero() const override;
  virtual bool has_inverse() const override;
  virtual bool is_cancellable() const override;

  virtual Polynomial<R> & operator + (const Ring_Element &) const override;
  virtual Polynomial<R> & operator - (const Ring_Element &) const override;
  virtual Polynomial<R> & operator * (const Ring_Element &) const override;

};

template<typename R>
Polynomial<R>::Polynomial() {
  deg = 0;
  coeffs.resize(1, 0);
}

template<typename R>
Polynomial<R>::Polynomial(const Polynomial<R> & other)
    : deg(other.deg), coeffs(deg + 1)
{
  for (unsigned i = 0; i <= deg; ++i)
    coeffs[i] = other.coeffs[i];
}

template<typename R>
Polynomial<R>::Polynomial(const initializer_list<R> & list)
    : deg(list.size() - 1), coeffs(deg + 1)
{
  DEGREE_TYPE i = 0;
  for (auto c : list) {
    coeffs[i] = c;
    ++i;
  }
  verify_degree();
}

template<typename R>
Polynomial<R>::Polynomial(const vector<R> & list)
    : deg(list.size()), coeffs(deg + 1)
{
  DEGREE_TYPE i = 0;
  for (auto c : list) {
    coeffs[i] = c;
    ++i;
  }
  verify_degree();
}

template<typename R>
void Polynomial<R>::verify_degree() {
  unsigned i;
  for (i = deg; i > 0 and coeffs[i].is_zero(); --i) { /* do nothing */ }
  deg = i;
}

template<typename R>
void Polynomial<R>::set_degree(DEGREE_TYPE d) {
  if (deg < d)
    coeffs.resize(d + 1);
  deg = d;
}

template<typename R>
void Polynomial<R>::clear_coeffs() {
  R zero;
  for (DEGREE_TYPE i = 0; i <= deg; ++i)
    coeffs[i] = zero;
}

template<typename R>
DEGREE_TYPE Polynomial<R>::degree() const { return deg; }

template<typename R>
const R & Polynomial<R>::coeff(DEGREE_TYPE d) const { return coeffs[d]; }

template<typename R>
const R & Polynomial<R>::operator() (const R & a) const {
  static R result;
  for (DEGREE_TYPE i = 0; i <= deg; ++i) {
    if (not coeffs[i].is_zero()) {
      R term(a);
      for (DEGREE_TYPE j = 0; j <= deg; ++j)
        term = term * a;
      term = term * coeffs[i];
      result = result + term;
    }
  }
  return result;
}

template<typename R>
void Polynomial<R>::set_coeff(DEGREE_TYPE d, const R & value) {
  coeffs[d] = value;
}

template<typename R>
const Polynomial<R> & Polynomial<R>::operator = (const Ring_Element & o) {
  if (*this != o) {
    const Polynomial<R> & other = dynamic_cast<const Polynomial<R> &>(o);
    coeffs.resize(other.degree() + 1);
    deg = other.degree();
    for (DEGREE_TYPE i = 0; i < deg; ++i)
      coeffs[i] = other.coeffs[i];
  }
  return *this;
}

template<typename R>
bool Polynomial<R>::operator == (const Ring_Element & o) const {
  auto & other = dynamic_cast<const Polynomial<R> &>(o);
  bool result = deg == other.deg;
  for (DEGREE_TYPE i = 0; result and i < deg; ++i)
    result = coeffs[i] == other.coeffs[i];
  return result;
}

template<typename R>
bool Polynomial<R>::operator != (const Ring_Element & o) const {
  auto & other = dynamic_cast<const Polynomial<R> &>(o);
  bool result = deg == other.deg;
  for (DEGREE_TYPE i = 0; result and i < deg; ++i)
    result = coeffs[i] == other.coeffs[i];
  return result;
}

template<typename R>
bool Polynomial<R>::operator == (const R & value) const {
  return deg == 0 and coeffs[0] == value;
}

template<typename R>
bool Polynomial<R>::is_one() const {
  return degree() == 0 and coeffs[0].is_one();
}

template<typename R>
bool Polynomial<R>::is_zero() const {
  return degree() == 0 and coeffs[0].is_zero();
}

template<typename R>
bool Polynomial<R>::has_inverse() const {
  return degree() == 0 and coeffs[0].is_cancellable();
}

template<typename R>
bool Polynomial<R>::is_cancellable() const {
  bool result = false;
  for (DEGREE_TYPE i = 0; i < deg and not result; ++i)
    result = not coeffs[i].is_zero() and coeffs[i].is_cancellable();
  return result;
}

template<typename R>
Polynomial<R> & Polynomial<R>::operator + (const Ring_Element & o)
    const
{
  static Polynomial<R> result;
  const Polynomial<R> & other = dynamic_cast<const Polynomial<R> &>(o);
  DEGREE_TYPE d = degree();
  if (d < other.degree()) d = other.degree();
  result.set_degree(d);
  result.clear_coeffs();
  DEGREE_TYPE i;
  for (i = 0; i <= degree() and i <= other.degree(); ++i)
    result.coeffs[i] = coeffs[i] + other.coeffs[i];
  for ( /* */ ; i <= degree(); ++i)
    result.coeffs[i] = coeffs[i];
  for ( /* */ ; i <= other.degree(); ++i)
    result.coeffs[i] = other.coeffs[i];
  result.verify_degree();
  return result;
}

template<typename R>
Polynomial<R> & Polynomial<R>::operator - (const Ring_Element & o)
    const
{
  static Polynomial<R> result;
  const Polynomial<R> & other = dynamic_cast<const Polynomial<R> &>(o);
  DEGREE_TYPE d = degree();
  if (d < other.degree()) d = other.degree();
  result.set_degree(d);
  result.clear_coeffs();
  DEGREE_TYPE i;
  for (i = 0; i <= degree() and i <= other.degree(); ++i)
    result.coeffs[i] = coeffs[i] - other.coeffs[i];
  for ( /* */ ; i <= degree(); ++i)
    result.coeffs[i] = coeffs[i];
  for ( /* */ ; i <= other.degree(); ++i)
    result.coeffs[i] = result.coeffs[i] - other.coeffs[i];
  result.verify_degree();
  return result;
}

template<typename R>
Polynomial<R> & Polynomial<R>::operator * (const Ring_Element & o)
    const
{
  static Polynomial<R> result;
  const Polynomial<R> & other = dynamic_cast<const Polynomial<R> &>(o);
  DEGREE_TYPE d = degree() + other.degree();
  result.set_degree(d);
  result.clear_coeffs();
  for (DEGREE_TYPE i = 0; i <= degree(); ++i)
    for (DEGREE_TYPE j = 0; j <= other.degree(); ++j)
      result.coeffs[i + j] = result.coeffs[i + j] + coeffs[i] * other.coeffs[j];
  result.verify_degree();
  return result;
}

template<typename R>
void Polynomial<R>::common_division(const Polynomial<R> & other) const {
  last_divisor = other;
  if (deg < other.deg) last_quotient.set_degree(0);
  else last_quotient.set_degree(deg - other.deg);
  last_quotient.clear_coeffs();
  last_remainder = *this;
  while (last_remainder.deg >= other.deg) {
    DEGREE_TYPE d = last_remainder.deg - other.deg;
    last_quotient.coeffs[d] = last_remainder.coeffs[last_remainder.deg];
    for (unsigned i = 0; i <= other.deg; ++i)
      last_remainder.coeffs[i + d] = last_remainder.coeffs[i + d]
          - other.coeffs[i]*last_remainder.coeffs[last_remainder.deg];
    last_remainder.verify_degree();
  }
}

template<typename R>
const Polynomial<R> & Polynomial<R>::operator / (const Polynomial<R> & other)
    const
{
  static Polynomial<R> result;
  if (other == last_divisor)
    result = last_quotient;
  else {
    common_division(other);
    result = last_quotient;
  }
  return result;
}

template<typename R>
const Polynomial<R> & Polynomial<R>::operator % (const Polynomial<R> & other)
    const
{
  static Polynomial<R> result;
  if (other == last_divisor)
    result = last_remainder;
  else {
    common_division(other);
    result = last_remainder;
  }
  return result;
}

template<typename R>
ostream & operator << (ostream & os, const Polynomial<R> & p) {
  bool first = true;
  for (DEGREE_TYPE d = p.degree(); d > 0; --d) {
    if (not p.coeff(d).is_zero()) {
      if (first) first = false;
      else {
        os << " + ";
      }
      if (not p.coeff(d).is_one()) os << p.coeff(d) << " ";
      os << "x";
      if (d > 1)
        os << "^" << d;
    }
  }
  if (p.degree() == 0 or not p.coeff(0).is_zero()) {
    if (not first) os << " + ";
    os << p.coeff(0);
  }
  return os;
}

template <typename R>
Polynomial<R> Polynomial<R>::last_divisor;
template <typename R>
Polynomial<R> Polynomial<R>::last_quotient;
template <typename R>
Polynomial<R> Polynomial<R>::last_remainder;

}

#endif