polynomial.h

This library defines basic operation on polynomials, and contains also 3 different roots (zeroes)-fi

#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
/** 
*       This program has been developed by Micha雔 Roy.  
*	You have the author's permission to incorporate this code 
*       into your product, royalty free.  

*       The author specifically disclaims all warranties, express or
*       implied, and all liability, including consequential and
*       other indirect damages, for the use of this code, 
*       including liability for infringement of any proprietary
*       rights, and including the warranties of merchantability
*       and fitness for a particular purpose.  And he does not 
*       assume any responsibility for any errors which may 
*       appear in this code nor any responsibility to update it.
*
*   (c) 2005, Michael Roy.  [e-mail]
*/
/** 
*   @file polynomial.h
*
*   Polynomial class definition
*
*   @ingroup Poly
*   @author Michael Roy <a href="mailto:[e-mail]">[e-mail]</a>
*   @version 1.0
*
*   @see    polyzero.h\n
*           laguerre.h\n
*           muller.h
*/

#include <vector>
#include <algorithm>
#include <complex>
#include <stdarg.h>
#include "floatspecs.h"

/** @defgroup Poly The Polynomial Templates
*   The Polynomial Templates.
*/
/** @defgroup PolyEval Polynomial evaluation functions
*   Polynomial evaluation functions.
*   @ingroup Poly
*/
/** @defgroup PolyOps Polynomial arithmetics functions
*   Polynomial arithmetics functions.
*   @ingroup Poly
*/
/** @addtogroup Poly 
*   @{ */
/**
*   defines operations on polynomials.
*
*   Class Polynomial is based on the class \p std::vector, and thus 
*   requires the Standard C++ Library (STL).
*
*   Coeficients are stored with lowest degree coeficent first, i.e.: the 
*   polynomial \f${a_0}+{a_1}x+{a_2}x^2+{a_3}x^3\f$ will be stored 
*   in the vector as \f$[a_0, a_1, a_2, a_3]\f$.
*
*   The coeficient \f$a_0\f$ of degree 0 is always populated.  That means
*   \p P[0] will always be valid for a polynomial \p P.
*
*   Evaluations are made using Horner's recurence, or synthetic division.
*
*   Valid types for template parameter \e T are any floating point,
*   or complex types which define operators +, -, *, /, and functions abs().
*
*   Several constructors are available, the following example shows them and 
*   their usage.
*
*   @par Example:
*   @include constructor.cpp
*
*   @author Michael Roy <a href="mailto:[e-mail]">[e-mail]</a>
*   @version 0.1
*
*   @see    polynomial.h for general support function\n
*           polyzero.h for polynomial roots support funtions\n
*           laguerre.h for roots-finding by Laguerre's method\n
*           muller.h for roots-finding by Muller's method
*/

template <class T>
class Polynomial : public std::vector<T> {

public:
    /** @name Construction */
    /** @{ */
    /**
    *   default constructor.
    *
    *   Creates an empty polynomial of degree zero, you can optionally specify
    *   a scalar value.
    *
    *   @param[in]  a   default scalar value for coeficient \f$a_0\f$.
    */

    explicit Polynomial(const T& a = T(0)) {
        push_back(a);
    }

    /**
    *   copy constructor.
    *
    *   Creates an exact copy of Polynomial \f$p\f$
    *
    *   @param[in]  p   source polynomial.
    *
    */

    explicit Polynomial(const std::vector<T>& p) : std::vector(p) { }

    /**
    *   initialization constructor.
    * 
    *   Initializes a polynomial of arbitrary size and coeficients.  Care must be taken
    *   to properly cast the variable arguments to the same type as template type \e T.
    *
    *   @param[in]  n       polynomial's degree
    *   @param[in]  a0      degree 0 coeficient
    *   @param[in]  ...     coeficients \f$a_1..a_n\f$, number of parameters must match 
    *                       parameter \a n.
    */

    Polynomial(unsigned n, T a0, ...) {
        va_list marker;
        int i = n;
        if (sizeof(T) != sizeof(float)) {
            push_back(a0);
            va_start(marker, a0);
            while (i-- > 0) {
                T an = va_arg(marker, T);
                push_back(an);
            }
            va_end(marker);
            compact();
        }
        else {
            // in c/c++, floats are pushed as doubles
            push_back(a0);
            va_start(marker, a0);
            while (i-- > 0) {
                double an = va_arg(marker, double);
                push_back(T(an));
            }
            va_end(marker);
            compact();
        }
    }
    /** @} construction */

    /** @name Methods */
    /** @{ */
    /**
    *   compacts array.
    *
    *   Compacts array size to match actual polynomial degree by enforcing 
    *   the rule <i>p[p.size() - 1] != 0</i>.
    *
    *   This function is called by other member functions to enforce the rule
    *   automatically whenever operations are performed on the polynomial.
    *
    *   @see degree()
    */

    void compact() {
        while((back() == T(0)) && size() > 1) pop_back();
    }
    
    /**
    *   returns degree of polynomial.
    *
    *   Degree of the polynomial is the index of the highest degree non_null
    *   coeficient.
    *
    *   @par Example:
    *   @include degree.cpp
    *   @par Program output:
    *   @include degree.txt
    *
    *   @returns the degree of the polynomial.
    *
    *   @see degree()
    */

    int degree() const {
        int m = (int)size() - 1;
        while ((m > 0) &&( (*this)[m] == T(0))) --m;
        return m;
    }

    /**
    *   returns degree of polynomial.
    *
    *   Degree of the polynomial is the index of the highest degree non_null
    *   coeficient.
    *
    *   This non-constant version of the function calls compact() before 
    *   returning.
    *
    *   @returns The degree of the polynomial.
    *
    *   @see degree() const
    */

    int degree() {
        compact();
        int m = (int)size() - 1;
        if (m < 0) {
            resize(1);
            m = 0;
        }
        return (int)size() - 1;
    }

    /**
    *   shifts the coeficients.
    *
    *   Shifts the polynomial coeficients left or right, effectively 
    *   doing the operation \f$q(x)=x^{n}p(x)\f$.
    *
    *   If \a n is negative, the remainder of the division 
    *   \f$\frac{p(x)}{x^{-n}}\f$ is lost.
    *
    *   @param[in]  n   power of multipler \f$x^n\f$.
    */

    void shift(int n) {
        if (n < 0) {
            int newSize = (int)size() + n;
            if (newSize >= 1) {
                std::copy(&(*this)[-n], &*end(), &(*this)[0]);
            }
            else {
                newSize = 1;
                (*this)[0] = T(0);
            }
            resize(newSize);
        }
        else if (n > 0) {
            resize(size() + n);
            std::copy(begin(), end(), begin() + n);
        }
    }

    /**
    *   makes polynomial monic.
    *
    *   Makes polynomial monic. A polynomial is said monic if its highest 
    *   degree coeficient equals 1.
    *
    *   @see isMonic() const
    */

    void makeMonic() {
        compact();
        if (!isMonic())
            *this /= (*this)[degree()];
    }

    /**
    *   indicates if polynomial is monic.
    *   Indicates if polynomial is monic. A polynomial is said monic if its 
    *   highest degree coeficient equals 1.
    *
    *   @returns \p true if polynomial is monic.
    *
    *   @see makeMonic()
    */

    inline bool isMonic() const {
        return ((*this)[degree()] == T(1));
    }
 
    /**
    *   get derivative of \f$p(x)\f$.
    *
    *  Computes the first derivative polynomial \f$p'(x)\f$ of polynomial \f$p(x)\f$. 
    *
    *   @param[out] pp    on exit, contains the first derivative \f$p'(x)\f$.
    *
    *   @see evalAndDerive(const T& x, T& ppx) const
    */

    void getDerivative(Polynomial<T>& pp) const {
        unsigned i;
        T ii;
        pp.resize(size() - 1);
        for (i = 0, ii = T(1); i < pp.size(); ++i, ii += T(1))
            pp[i] = ii * (*this)[i + 1];
    }

    /** @} methods */

    /** @name Evaluation */
    /** @{ */
    /**
    *   evaluates polynomial value \f$p(x)\f$.
    *
    *   Evaluates \f$p(x)\f$ by Horner's recurence.
    *
    *   No error checking is performed, so any validation for overflow, 
    *   or underflow is the responsibility of the caller.
    *
    *   @param[in]  x   parameter of \f$p(x)\f$
    *
    *   @returns The value of \f$p(x)\f$.
    */

    T eval(const T& x) const {
        return ::eval(*this, x);
    }

    /**
    *   evaluates polynomial value \e p(a), and computes \f$\frac{p(x)}{(x - a)}\f$.
    *
    *   Evaluates \f$p(a)\f$ and polynomial \f$q(x)=\frac{p(x)-p(a)}{(x - a)}\f$ at the same time 
    *   by Horner's recurence.
    *
    *   It can be shown that \f$p(x)=q(x)(x - a)+p(a)\f$, 
    *   which essentially means that \f$p(a)\f$ is the remainder of 
    *   \f$\frac{p(x)}{(x - a)}\f$. For the special case where \f$a\f$ is a root of 
    *   \f$p(x)\f$ (i.e. \f$p(a)=0\f$), this operation is called 
    *   <i>deflation</i>.  That is, on exit \f$q(x)\f$ will have all the roots of 
    *   \f$p(x)\f$, except for the root \f$a\f$.
    *
    *   No error checking is performed, so any validation for overflow, 
    *   or underflow is the responsibility of the caller.
    *
    *   @param[in]  a   parameter of \f$p(a)\f$.
    *   @param[out] q   on exit contains the quotient \f$q(x)=\frac{p(x)-p(a)}{(x - a)}\f$.
    *   @returns \f$p(a)\f$, the remainder of \f$\frac{p(x)}{(x - a)}\f$.
    *
    *   @see evalAndDeflate(const Polynomial<T>& p, const U& x, Polynomial<U>& q, V& e)
    */

    T evalAndDeflate(const T& a, Polynomial<T>& q) const {
        return ::evalAndDeflate(*this, a, q);
    }

    /**
    *   evaluates \f$p(x)\f$ and derivate \f$dy=p'(x)\f$ at the same time.
    *
    *   Evaluates efficiently \f$p(x)\f$ and first derivative \f$p'(x)\f$
    *   at the same time.
    *
    *   @param[in]  x   variable to evaluate.
    *   @param[out] ppx on exit contains \f$p'(x)\f$.
    *   @returns The value of \f$p(x)\f$.
    */

    T evalAndDerive(const T& x, T& ppx) const {
        int m = degree();
        T y;
        ppx = 0;
        y = (*this)[m];
        for (int i = m - 1; i >= 0; --i) {
            ppx = (ppx * x) + y;
            y = (y * x) + (*this)[i];
        }
        return y;
    }
    /** @} evaluation */

    /** @name Operators */
    /** @{ */
    /**
    *   make a copy of a polynomial.
    *
    *   Make a copy of a polynomial.
    *
    *   @param[in]  q   polynomial to copy.
    *
    *   @returns A const reference on this object.
    */

    inline const Polynomial<T>& operator= (const std::vector<T>& q) {
        (std::vector<T>&)*this = q;
        return *this;
    }

    /**
    *   in-place polynomial addition.
    *
    *   Performs the addition \f$r=p+q\f$, storing the result stored in \p this.
    *
    *   @param[in]  q   a polynomial.
    *   @returns        A const reference to \p this.
    *
    *   @see    operator+ (const Polynomial<T>& p, const Polynomial<T>& q)\n
    *           Polynomial::operator-= (const Polynomial<T>& A)\n
    *           add(Polynomial<T>& p, const Polynomial<T>& q)
    */

    inline const Polynomial<T>& operator+= (const Polynomial<T>& q) {
        add(*this, q);
        return *this;
    }

    /**
    *   in-place polynomial substraction.
    *   Performs the substraction \f$r=p-q\f$ and stores the resulting polynomial in 
    *   \p this.
    *
    *   @param[in]  q   a polynomial.
    *   @returns        A const reference to \p this.
    *
    *   @see    operator- (const Polynomial<T>& p, const Polynomial<T>& q)\n
    *           Polynomial::operator+= (const Polynomial<T>& A)\n
    *           sub(Polynomial<T>& p, const Polynomial<T>& q)
    */

    inline const Polynomial<T>& operator-= (const Polynomial<T>& q) {
        sub(*this, q);
        return *this;
    }

    /**
    *   in-place polynomial multiply.
    *
    *   Multiplies two polynomials by convolution.  Performs the operation 
    *   \f$p\leftarrow pq\f$
    *
    *   @param[in]  p   a polynomial to multiply to \p this object.
    *   @returns        A const reference to this.
    *
    *   @see    mul(Polynomial<T> & r, const Polynomial<T>& p, const Polynomial<T>& q)\n
    *           Polynomial::operator* (const Polynomial<T>& p, const Polynomial<T>& q)
    */

    inline const Polynomial<T> operator*= (const Polynomial<T>& p) {
        Polynomial<T> q = *this;
        mul(*this, p, q);
        return *this;
    }

    /**
    *   in-place polynomial scaling.
    *
    *   Polynomial scaling, performs the operation \f$p(x)\leftarrow s.p(x)\f$.