polynomial.h

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

    *
    *   @returns A constant reference to \p this.
    *
    *   @see    mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p)\n
    *           Polynomial<T>::operator*= (const T& f)\n
    *           operator/ (const T& f, const Polynomial<T>& p)
    */

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

    /**
    *   in-place polynomial scaling.
    *
    *   Polynomial scaling, performs the operation \f$p(x)\leftarrow \frac{p(x)}{s}\f$.
    *
    *   @returns A constant reference to \p this.
    *
    *   @see    mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p)\n
    *           Polynomial<T>::operator*= (const T& s)\n
    *           operator/ (const Polynomial<T>& p, const T& s)
    */

    Polynomial<T>& operator/= (const T& s) {
        mul(*this, T(1) / s);
        return *this;
    }
    /** @} operators */

    /** @name Non-member operators */
    /** @{ */
    /**
    *   out-of-place polynomial addition.
    *
    *   Performs the addition \f$r=p+q\f$.
    *
    *   Calling directly the three parameters function add() is actually more 
    *   efficient, as it will avoid an unnecessary copy of the polynomial.
    *
    *   @see    operator- (const Polynomial<T>& p, const Polynomial<T>& q)\n
    *           Polynomial<T>::operator+= (const Polynomial<T>& p)\n
    *           add(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q)
    */

    friend inline Polynomial<T> operator+ (const Polynomial<T>& p, const Polynomial<T>& q) {
        Polynomial<T> r;
        add(r, p, q);
        return r;
    }

    /**
    *   out-of-place polynomial substraction.
    *
    *   Performs the substraction \f$r=p-q\f$.
    *
    *   @param[in]  p   a polynomial.
    *   @param[in]  q   another polynomial.
    *   @returns    The result @e r of the operation \f$r=p-q\f$.
    *
    *   @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)
    */

    friend inline Polynomial<T> operator- (const Polynomial<T>& p, const Polynomial<T>& q) {
        Polynomial<T> r;
        sub(r, a, b);
        return r;
    }

    /**
    *   out-of-place polynomial multiply.
    *
    *   Multiplies two polynomials by convolution.  It is advisable to use 
    *   mul() instead, as it avoids the copying of the resulting polynomial
    *   off the stack.
    *
    *   @param[in]  p   a polynomial.
    *   @param[in]  q   another polynomial.
    *   @returns        The product \f$r=pq\f$
    *
    *   @see    mul(Polynomial<T> & r, const Polynomial<T>& p, const Polynomial<T>& q)\n
    *           Polynomial::operator*= (const Polynomial<T>& p)
    */

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

    /**
    *   out-of_place polynomial scaling.
    *
    *   Polynomial scaling.
    *
    *   @returns The polynomial \f$r(x)=s.p(x)\f$.
    *
    *   @see    mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p)\n
    *           Polynomial<T>::operator*= (const T& s)\n
    *           operator/ (const T& s, const Polynomial<T>& p)
    */

    friend inline Polynomial<T> operator* (const T& s, const Polynomial<T>& p) {
        Polynomial<T> r;
        mul(r, s, p);
        return r;
    }

    /**
    *   out-of-place polynomial scaling.
    *
    *   Polynomial scaling, performs the operation \f$r(x)=\frac{p(x)}{s}\f$.
    *
    *   @returns A constant reference to \p this.
    *
    *   @see    mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p)\n
    *           Polynomial<T>::operator*= (const T& s)\n
    *           operator/= (const T& s)
    */

    friend Polynomial<T> operator/ (const Polynomial<T>& p, const T& s) {
        Polynomial<T> r;
        mul(r, T(1) / f, p);
        return r;
    }
    /** @} Non-member operators */
};

/**
*   out-of-place polynomial addition.
*
*   Performs addition of two polynomials \f$r=p+q\f$, 
*   resulting polynomial is stored in user-supplied polynomial \a r.
*
*   @param[out] r   on exit, contains the polynomial \f$r=p+q\f$.
*   @param[in]  p   first input polynomial.
*   @param[in]  q   second input polynomial.
*
*   @ingroup PolyOps
*   @see    operator+ (const Polynomial<T>& p, const Polynomial<T>& q)\n
*           add(Polynomial<T>& p, const Polynomial<T>& q)
*/

template <class T>
void add(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q) {
    int i, m, n;
    m = p.degree();
    n = q.degree();
    if (m <= n) {
        r.resize(n + 1);
        for (i = 0; i <= m; ++i)
            r[i] = p[i] + q[i];
        for ( ; i <= n; ++i)
            r[i] = q[i];
    }
    else {
        r.resize(m + 1);
        for (i = 0; i <= n; ++i)
            r[i] = p[i] + q[i];
        for ( ; i <= m; ++i)
            r[i] = p[i];
    }
    r.compact();
}

/**
*   in-place polynomial addition.
*
*   Performs in-place addition of two polynomials \f$p\leftarrow p+q\f$.
*
*   @param[in,out]  p   first input polynomial.  On exit, contains the sum \f$p+q\f$.
*   @param[in]      q   second input polynomial.
*
*   @ingroup PolyOps
*   @see    add(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q)\n
*           Polynomial::operator += (const Polynomial<T>& A)
*/

template <class T>
void add(Polynomial<T>& p, const Polynomial<T>& q) {
    int m = p.degree();
    int n = q.degree();
    if (m < n)
        p.resize(n + 1, T(0));

    for (int i = 0; i <= n; ++i)
        p[i] += q[i];
    p.compact();
}

/**
*   out-of-place substraction of polynomials.
*
*   Performs substraction of two polynomials \f$r=p-q\f$, 
*   resulting polynomial is stored in user-supplied polynomial \a r.
*
*   @param[out] r   on exit, contains the polynomial \f$p-q\f$.
*   @param[in]  p   first input polynomial.
*   @param[in]  q   second input polynomial.
*
*   @ingroup PolyOps
*/

template <class T>
void sub(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q) {
    int i;
    int m = p.degree();
    int n = q.degree();
    if (m <= n) {
        r.resize(n + 1);
        for (i = 0; i <= m; ++i)
            r[i] = p[i] - q[i];
        for ( ; i <= n; ++i)
            r[i] = -q[i];
    }
    else {
        r.resize(m + 1);
        for (i = 0; i <= n; ++i)
            r[i] = p[i] - q[i];
        for ( ; i <= m; ++i)
            r[i] = p[i];
    }
    r.compact();
}

/**
*   in-place substraction of polynomials.
*
*   Performs in-place substraction of two polynomials \f$p\leftarrow p-q\f$.
*
*   @param[in,out]  p   first input polynomial.  On exit, contains the difference \f$p-q\f$.
*   @param[in]      q   second input polynomial.
*
*   @ingroup PolyOps
*/

template <class T>
void sub(Polynomial<T>& p, const Polynomial<T>& q) {
    int i;
    int m = p.degree();
    int n = q.degree();
    if (m < n)
        p.resize(n + 1, T(0));
    for (i = 0; i <= n; ++i)
        p[i] -= q[i];
    p.compact();
}

/**
*   out-of-place polynomial division.
*
*   Performs polynomial division.  On exit:
*   \f[q=\frac{u}{v}\f]
*   \f[r=u\%v\f]
*
*   @param[out] q   quotient \f$q=\frac{u}{v}\f$
*   @param[out] r   remainder \f$r=u\%v\f$
*   @param[in]  u   numerator.
*   @param[in]  v   divisor.
*
*   @ingroup PolyOps
*   @see div(Polynomial<T>& u, Polynomial<T>& r, const Polynomial<T>& v)
*/

template <class T>
void div(Polynomial<T>& q, Polynomial<T>& r, const Polynomial<T>& u, const Polynomial<T>& v) {
    int na = u.degree();
    int nb = v.degree();
    if (na >= nb) {
        r = u;
        T t;
        q.resize(1 + na - nb);
        std::fill(&q[0], &q[q.size()], T(0));
        for (int i = na - nb; i >= 0; --i) {
            q[i] = r[nb + i] / v[nb];
            r[nb + i] = T(0);
            for (int j = nb + i - 1; j >= i; --j) {
                r[j + 1] -= (q[i] * v[j - i]);
            }
            r[i + nb] = 0;
        }
        r.compact();
    }
    else {
        r = u;
        q.resize(1);
        q[0] = T(0);
    }
}

/**
*   in-place polynomial division.
*
*   Performs polynomial division.  On exit:
*   \f[u\leftarrow \frac{u}{v}\f]
*   \f[r=u\%v\f]
*
*   @param[in,out]  u   numerator, holds quotient \f$\frac{u}{v}\f$ on exit.
*   @param[out]     r   remainder \f$r=u\%v\f$
*   @param[in]      v   divisor.
*
*   @ingroup PolyOps
*   @see div(Polynomial<T>& q, Polynomial<T>& r, const Polynomial<T>& u, const Polynomial<T>& v)
*/

template <class T>
void div(Polynomial<T>& u, Polynomial<T>& r, const Polynomial<T>& v) {
    int na = u.degree();
    int nb = v.degree();
    if (na >= nb) {
        r = u;
        u.resize(na - nb + 1);
        for (int i = na - nb; i >= 0; --i) {
            u[i] = r[nb + i] / v[nb];
            r[nb + i] = T(0);
            for (int j = nb + i - 1; j >= i; --j)
                r[j] -= (u[i] * v[j - i]);
        }
        r.compact();
    }
    else {
        r = u;
        u.resize(1);
        u[0] = T(0);
    }
}

/**
*   out-of-place polynomial scaling.
*
*   Polynomial scaling. Performs the operation \f$r(x)=s.p(x)\f$.
*
*   @param[out] r   scaled polynomial \f$r(x)=s.p(x)\f$
*   @param[in]  s   scalar.
*   @param[in]  p   A polynomial.
*
*   @ingroup PolyOps
*   @see    mul(Polynomial<T>& p, const T& s)\n
*           operator* (const T& s, const Polynomial<T>& p)\n
*           operator* (const Polynomial<T>& p, const T& s)
*/

template <class T>
void mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p) {
    if (s != T(1)) {
        r.resize(p.degree() + 1);
        for (unsigned i = 0; i < r.size(); ++i)
            r[i] = s * p[i];
    }
}

/**
*   in-place polynomial scaling.
*
*   Polynomial scaling. Performs the operation \f$p(x)\leftarrow s.p(x)\f$.
*
*   @param[in, out] p   A polynomial, on exit contains the product \f$s.p(x)\f$
*   @param[in]      s   scalar.
*
*   @ingroup PolyOps
*   @see    mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p)\n
*           Polynomial<T>::operator*= (Polynomial<T>& p, const T& s)
*/

template <class T>
inline void mul(Polynomial<T>& p, const T& s) {
    mul(p, s, p);
}

/**
*   out-of-place polynomial multiply.
*
*   Multiplies two polynomials by convolution.
*
*   @param[out] r   resulting polynomial \f$r=pq\f$
*   @param[in]  p   a polynomial.
*   @param[in]  q   another polynomial.
*
*   @ingroup PolyOps
*   @see    operator* (const Polynomial<T>& p, const Polynomial<T>& q)\n
*           Polynomial<T>::operator*= (const Polynomial<T>& p)
*/

template <class T>
void mul(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q) {
    int na = p.degree();
    int nb = q.degree();
    r.resize(na + nb + 1);
    std::fill(&r[0], &r[r.size()], T(0));
    for (int i = 0; i <= na; ++i)
        for (int j = 0; j <= nb; ++j)
            r[i + j] += p[i] * q[j];
}

/**
*   evaluates polynomial value \f$y = p(x)\f$.
*
*   Evaluates \f$p(x)\f$ by Horner's recurence.  This is the preferred way of 
*   evaluating polynomials, since the member function Polynomial<T>::eval() is
*   resricted to same-type polynomial coeficients and variable \a x
*
*   No error checking is performed, so any validation for overflow, 
*   or underflow is the responsibility of the caller.
*
*   @param[in]  p   a polynomial.
*   @param[in]  x   parameter of \f$p(x)\f$
*   @returns The value of \f$p(x)\f$.
*
*   @ingroup PolyEval
*   @see    Polynomial<T>::eval(const T& x) const\n
*           eval(const Polynomial<T>& p, const U& x, V& e)
*/

template<class T, class U>
U eval(const Polynomial<T>& p, const U& x) {
    if (x == U(0)) 
        return p[0];
    else {
        int i, m = p.degree();
        register U px = p[m];
        for (i = m - 1; i >= 0; --i) {
            px *= x;
            px += p[i];
        }
        return px;
    }
}

/**
*   evaluates polynomial value \f$y = p(x)\f$.
*
*   Evaluates \f$p(x)\f$ by Horner's recurence.  This is the preferred way of 
*   evaluating polynomials, since the member function Polynomial<T>::eval() is
*   resricted to same-type polynomial coeficients and variable \a x
*
*   No error checking is performed, so any validation for overflow, 
*   or underflow is the responsibility of the caller.
*
*   @param[in]  p   a polynomial.