Polynomial.h

#ifndef UTIL_POLYNOMIAL_H
#define UTIL_POLYNOMIAL_H
 
/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <util/containers/GArray.h>  // base class
#include <util/math/Rational.h>      // default template argument
#include <util/math/Binomial.h>      // used in implementation
#include <util/containers/DArray.h>  // used in implementation
#include <util/global.h>
 
 
#include <iostream>
 
namespace Util
{
 
   template <typename T = Rational>
   class Polynomial : private GArray<T>
   {
 
   public:
 
 
 
      Polynomial(int capacity = 10);
 
      explicit Polynomial(T c);
 
      Polynomial(Array<T> const & coeffs);
 
      Polynomial(Polynomial<T> const & other);
 
      template <typename U>
      Polynomial<T>& operator = (Polynomial<U> const & other);
 
      void setToZero();
 
 
 
      // Get a specific coefficient by index.
      using GArray<T>::operator [];
 
      int degree() const;
 
      // Logical size of the coefficient array.
      using GArray<T>::size;
 
      // Physical capacity of the coefficient array.
      using GArray<T>::capacity;
 
 
 
      Polynomial<T>& operator += (Polynomial<T> const & a);
 
      Polynomial<T>& operator += (T a);
 
      Polynomial<T>& operator -= (Polynomial<T> const & a);
 
      Polynomial<T>& operator -= (T a);
 
      Polynomial<T>& operator *= (T a);
 
      Polynomial<T>& operator /= (T a);
 
      Polynomial<T>& operator *= (Polynomial<T> const & a);
 
 
 
      Polynomial<T> integrate() const;
 
      Polynomial<T> differentiate() const;
 
      Polynomial<T> reflect() const;
 
      Polynomial<T> shift(T a) const;
 
 
 
      T operator () (T x) const;
 
      double evaluate (double x) const;
 
 
      // Static public member functions
      
      static Polynomial<T> monomial(int n);
 
   };
 
   // Inline methods
 
   /*
   * Constructor.
   */
   template <typename T>
   inline 
   Polynomial<T>::Polynomial(int capacity)
   {  GArray<T>::reserve(capacity); }
 
   /*
   * Construct a constant polynomial f(x) = c.
   */
   template <typename T>
   inline
   Polynomial<T>::Polynomial(T c)
   {
      GArray<T>::reserve(10); 
      GArray<T>::append(c);
   }
 
   /*
   * Construct from array of coefficients.
   */
   template <typename T>
   inline
   Polynomial<T>::Polynomial(Array<T> const & coeffs)
   {
      if (coeffs.capacity() > 0) {
         GArray<T>::reserve(coeffs.capacity());
         for (int i = 0; i < coeffs.capacity(); ++i) {
            GArray<T>::append(coeffs[i]);
         }
      }
   }
 
   /*
   * Copy constructor
   */
   template <typename T>
   inline
   Polynomial<T>::Polynomial(Polynomial<T> const & other)
   {
      GArray<T>::reserve(other.capacity());
      if (other.size() > 0) {
        for (int i = 0; i < other.size(); ++i) {
            GArray<T>::append(other[i]);
         }
      }
   }
 
   /*
   * Assignment from another polynomial.
   */
   template <typename T>
   template <typename U>
   inline
   Polynomial<T>& Polynomial<T>::operator = (Polynomial<U> const & other)
   {
      setToZero();
      if (other.size() >= 0) {
         if (other.size() > capacity()) {
            GArray<T>::reserve(other.capacity());
         }
         T element;
         for (int i = 0; i < other.size(); ++i) {
            element = other[i];
            GArray<T>::append(element);
         }
      }
      return *this;
   }
 
   /*
   * Set this polynomial to zero.
   */
   template <typename T>
   inline 
   void Polynomial<T>::setToZero()
   {  GArray<T>::clear(); }
 
   /*
   * Return degree of polynomial (size of coeff. array - 1).
   */
   template <typename T>
   inline
   int Polynomial<T>::degree() const
   {  return GArray<T>::size() - 1; }
 
   /*
   * Addition assignment operator : add another polynomial to this one.
   */
   template <typename T>
   Polynomial<T>& Polynomial<T>::operator += (Polynomial<T> const & a)
   {
      if (a.size() > 0) {
         int min = a.size() > size() ? size() : a.size();
         if (min > 0) {
            for (int i = 0; i < min; ++i) {
               (*this)[i] += a[i];
            }
         }
         if (a.size() > size()) {
            UTIL_CHECK(min == size());
            for (int i = size(); i < a.size(); ++i) {
               GArray<T>::append(a[i]);
            }
         }
      }
      return *this;
   }
 
   /*
   * Add a constant to this polynomial.
   */
   template <typename T>
   Polynomial<T>& Polynomial<T>::operator += (T a)
   {
      if (size() == 0) {
         (*this).append(a);
      } else {
         (*this)[0] += a;
      }
      return *this;
   }
 
   /*
   * Subtract assignment operator : subtract another polynomial from this.
   */
   template <typename T>
   Polynomial<T>& Polynomial<T>::operator -= (Polynomial<T> const & a)
   {
      if (a.size() > 0) {
         int min = a.size() > size() ? size() : a.size();
         if (min > 0) {
            for (int i = 0; i < min; ++i) {
               (*this)[i] -= a[i];
            }
         }
         if (a.size() > size()) {
            UTIL_CHECK(min == size());
            for (int i = size(); i < a.size(); ++i) {
               GArray<T>::append(-a[i]);
            }
         }
      }
      return *this;
   }
 
   /*
   * Subtract a constant from this polynomial.
   */
   template <typename T>
   Polynomial<T>& Polynomial<T>::operator -= (T a)
   {
      if (size() == 0) {
         (*this).append(-a);
      } else {
         (*this)[0] -= a;
      }
      return *this;
   }
 
   /*
   * Multipication assignment operator : multiply this by a scalar.
   */
   template <typename T>
   inline
   Polynomial<T>& Polynomial<T>::operator *= (T a)
   {
      if (size() >= 0) {
         for (int i = 0; i < size(); ++i) {
            (*this)[i] *= a;
         }
      }
      return *this;
   }
 
   /*
   * Division by scalar assignment operator.
   */
   template <typename T>
   inline
   Polynomial<T>& Polynomial<T>::operator /= (T a)
   {
      if (size() >= 0) {
         for (int i = 0; i < size(); ++i) {
            (*this)[i] /= a;
         }
      }
      return *this;
   }
 
   /*
   * Multiplication assignment operator : multiply this by a polynomial.
   */
   template <typename T>
   Polynomial<T>& Polynomial<T>::operator *= (Polynomial<T> const & a)
   {
      // If this polynomial is zero (size == 0), do nothing. Otherwise:
      if (size() > 0) {
 
         if (a.size() == 0) {
            // If polynomial a is zero, set this one to zero (clear it).
            setToZero();
         } else {
 
            // Compute size of new array of coefficients (degree + 1)
            int n = size() + a.size() - 1;
 
            // Make a copy of coefficients of this polynomial
            GArray<T> b(*this);
 
            // Clear this array of coefficients and reserve enough space
            setToZero();
            if (n > capacity()) {
               GArray<T>::reserve(n);
            }
 
            // Set all coefficients of resized array to zero
            int i;
            T zero = 0;
            for (i = 0; i < n; ++i) {
               GArray<T>::append(zero);
            }
 
            // Compute new coefficients as a double sum
            int j, k;
            for (i = 0; i < a.size(); ++i) {
               for (j = 0; j < b.size(); ++j) {
                  k = i + j;
                  UTIL_ASSERT(k < n);
                  (*this)[k] += a[i]*b[j];
               }
            }
 
         }
      }
      return *this;
   }
 
   /*
   * Compute and return indefinite integral of this polynomial.
   */
   template <typename T>
   Polynomial<T> Polynomial<T>::integrate() const
   {
      if (size() == 0) {
         Polynomial<T> b;
         b.setToZero();
         return b;
      } else {
         // Construct and compute coefficient array for integral
         DArray<T> coeffs;
         coeffs.allocate(size()+1);
         coeffs[0] = T(0);
         for (int i = 0; i < size(); ++i) {
            coeffs[i+1] = (*this)[i];
            coeffs[i+1] /= T(i+1);
         }
 
         // Construct and return associated Polynomial
         Polynomial<T> b(coeffs);
         return b;
      }
   }
 
   /*
   * Compute and return the derivatvie of this polynomial.
   */
   template <typename T>
   Polynomial<T> Polynomial<T>::differentiate() const
   {
      if (size() <= 1) {
 
         // If this polynomial is null or constant, return null
         Polynomial<T> b;
         b.setToZero();
         return b;
 
      } else {
 
         // Construct coefficient array for derivative polynomial
         DArray<T> coeffs;
         coeffs.allocate(size()-1);
         for (int i = 1; i < size(); ++i) {
            coeffs[i-1] = (*this)[i];
            coeffs[i-1] *= T(i);
         }
 
         // Construct and return associated Polynomial
         Polynomial<T> b(coeffs);
         return b;
      }
   }
 
   /*
   * Compute and return reflection g(x) = f(-x) this polynomial f(x).
   */
   template <typename T>
   Polynomial<T> Polynomial<T>::reflect() const
   {
      // Make copy
      Polynomial<T> b(*this);
 
      // Reverse odd power coefficients
      for (int i = 0; i < b.size(); ++i) {
         if (i%2 != 0) {
            b[i] *= -1;
         }
      }
 
      return b;
   }
 
   /*
   * Compute and return reflection g(x) = f(-x) this polynomial f(x).
   */
   template <typename T>
   Polynomial<T> Polynomial<T>::shift(T a) const
   {
      // Make copy of this
      Polynomial<T> b(*this);
 
      int degree = size() - 1;
      if (degree > 0) {
         Binomial::setup(degree);
         int n, m;
         T p;
         for (n = 1; n <= degree; ++n) {
            p = b[n]*a;
            for (m = 1; m <= n; ++m) {
               b[n -m] += Binomial::coeff(n, m)*p;
               p *= a;
            }
         }
      }
 
      return b;
   }
 
   /*
   * Evaluate polynomial at specific argument.
   */
   template <typename T>
   inline T Polynomial<T>::operator () (T x) const
   { 
      if (size() > 0) {
         int degree = size() - 1;
         T value = (*this)[degree];
         if (degree > 0) {
            for (int i = degree-1; i >= 0; --i) {
              value *= x;
              value += (*this)[i];
            }
         }
         return value;
      } else {
         return T(0);
      }
   }
 
   /*
   * Evaluate polynomial at specific floating point argument.
   */
   template <typename T>
   inline double Polynomial<T>::evaluate (double x) const
   {
      if (size() > 0) {
         int degree = size()-1;
         double value = (double)(*this)[degree];
         if (degree > 0) {
            for (int i = degree-1; i >= 0; --i) {
              value *= x;
              value += (double)(*this)[i];
            }
         }
         return value;
      } else {
         return 0.0;
      }
   }
 
   // Static function
 
   template <typename T>
   Polynomial<T> Polynomial<T>::monomial(int power)
   {
      UTIL_CHECK(power >= 0);
      Polynomial<T> a(power+1);
 
      // Append zero coefficients.
      T zero(0);
      for (int i = 0; i <= power; ++i) {
         a.GArray<T>::append(zero);
      }
 
      // Set coefficient of highest power to unity
      a[power] = T(1);
 
      return a;
   }
 
   // Related non-member functions
 
   template <typename T>
   bool operator == (Polynomial<T> const & a, Polynomial<T> const & b)
   {
      if (a.size() != b.size()) return false;
      if (a.size() > 0) {
         for (int i = 0; i < a.size(); ++i) {
            if (a[i] != b[i]) return false;
         }
      }
      return true;
   }
 
   template <typename T>
   bool operator != (Polynomial<T> const & a, Polynomial<T> const & b)
   {  return !(a == b); }
 
   template <typename T>
   inline
   Polynomial<T> operator - (Polynomial<T> const & a)
   {
      Polynomial<T> b(a);
      if (a.size() > 0) {
         for (int i = 0; i < a.size(); ++i) {
           b[i] *= -1;
         }
      }
      return b;
   }
 
   template <typename T>
   std::ostream& operator << (std::ostream& out, Polynomial<T> const & p)
   {
      out << "( ";
      if (p.size() > 0) {
         for (int i = 0; i < p.size(); ++i) {
            out << p[i] << " ";
         }
      }
      out << ")";
      return out;
   }
}
#endif