polynomial.h

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

*   @param[in]  x   parameter of \f$p(x)\f$
*   @param[out] e   round-off error to be expected in the calculation of 
*                   \f$p(x)\f$
*   @returns The value of \f$p(x)\f$.
*
*   @ingroup PolyEval
*   @see    eval(const Polynomial<T>& p, const U& x)\
*           Polynomial<T>::eval(const T& x) const
*/

template<class T, class U, class V>
U eval(const Polynomial<T>& p, const U& x, V& e) {
    static FloatSpecs<V> fpSpecs;
    if (x == U(0)) {
        e = V(0);
        return p[0];
    }
    else {
        int i, m = p.degree();
        if (m >= 1) {
            U px = p[m];
            V mx, mpx;
            mx = abs(x);
            e = (abs(px) * fpSpecs.MRE) / (fpSpecs.ARE + fpSpecs.MRE);
            for (i = m - 1; i >= 0; --i) {
                px *= x;
                px += p[i];
                //px = (px * x) + p[i];
                mpx = abs(px);
                e = (mx * e) + mpx;
            }
            e = (e * (fpSpecs.ARE + fpSpecs.MRE)) - (mpx * fpSpecs.MRE);
            return px;
        }
        else {
            e = V(0);
            return p[0];
        }
    }
}

/**
*   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]  p   a polynomial.
*   @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$.
*
*   @ingroup PolyEval
*   @see    Polynomial<T>::evalAndDeflate(const U& x, Polynomial<U>& q) const\n
*           Polynomial<T>::evalAndDeflate(const U& x, Polynomial<U>& q, V& e) const\n
*           evalAndDeflate(const Polynomial<T>& p, const U& x, Polynomial<U>& q, V& e)
*/

template<class T, class U>
U evalAndDeflate(const Polynomial<T>& p, const U& a, Polynomial<U>& q) {
    int i, m = p.degree();
    if (m > 0) {
        q.resize(m);
        if (a == T(0)) {
            std::copy(&p[1], &p[m + 1], &q[0]);
            return p[0];
        }
        else {
            U pa;
            pa = p[m];
            for (i = m - 1; i >= 0; --i) {
                q[i] = pa;
                pa *= a;
                pa += p[i];
                //pa = (pa * a) + p[i];
            }
            return pa;
        }
    }
    else {
        q.resize(1);
        q[0] = 0;
        return p[0];
    }
}

/**
*   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]  p   a polynomial.
*   @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$.
*   @param[out] e   round-off error to be expected in the calculation of 
*                   \f$p(a)\f$
*   @returns \f$p(a)\f$, the remainder of \f$\frac{p(x)}{(x - a)}\f$.
*
*   @ingroup PolyEval
*   @see    Polynomial<T>::evalAndDeflate(const U& x, Polynomial<U>& q) const\n
*           Polynomial<T>::evalAndDeflate(const U& x, Polynomial<U>& q, V& e) const\n
*           evalAndDeflate(const Polynomial<T>& p, const U& x, Polynomial<U>& q)
*/

template<class T, class U, class V>
U evalAndDeflate(const Polynomial<T>& p, const U& a, Polynomial<U>& q, V& e) {
    static FloatSpecs<V> fpSpecs;
    int i, m = p.degree();
    if (m > 0) {
        q.resize(m);
        if (a == T(0)) {
            e = 0;
            std::copy(&p[1], &p[m + 1], &q[0]);
            return p[0];
        }
        else {
            U pa;
            V mx, mpa;
            pa = p[m];
            mx = abs(a);
            e = (abs(pa) * fpSpecs.MRE) / (fpSpecs.ARE + fpSpecs.MRE);
            for (i = m - 1; i >= 0; --i) {
                q[i] = pa;
                pa *= a;
                pa += p[i];
                mpa = abs(pa);
                e = (mx * e) + mpa;
            }
            e = (e * (fpSpecs.ARE + fpSpecs.MRE)) - (mpa * fpSpecs.MRE);
            return pa;
        }
    }
    else {
        e = 0;
        q.resize(1);
        q[0] = 0;
        return p[0];
    }
}

/**
*   evaluates round-off error.
*
*   Evaluates round_off error in computing \f$p(x)\f$.
*
*   @param[in]  p   a polynomial.
*   @param[in]  mx  absolute value of x.
*   @returns    e   round-off error to be expected in the calculation of 
*                   \f$p(x)\f$
*
*   @ingroup PolyEval
*   @see    eval(const Polynomial<T>& p, const U& x, V& e)\n
*           evalAndDeflate(const Polynomial<T>& p, const U& x, Polynomial<U>& q, V& e)\n
*           Polynomial<T>::evalAndDeflate(const U& x, Polynomial<U>& q, V& e) const
*/

template <class T, class U>
U evalError(const Polynomial<T>& p, const U& mx) {
    static FloatSpecs<U> fpSpecs;
    T px;
    U e, mpx;
    int m = p.degree();
    px = p[m];
    e = (abs(px) * fpSpecs.MRE) / (fpSpecs.ARE + fpSpecs.MRE);
    for (int i = m - 1; i >= 0; --i) {
        px *= mx;
        px += p[i];
        mpx = abs(px);
        e = (mx * e) + mpx;
    }
    return (e * (fpSpecs.ARE + fpSpecs.MRE)) - (mpx * fpSpecs.MRE);
}

/**
*   evaluates \f$p(x)\f$ and first 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.
*
*   @ingroup PolyEval
*   @param[in]  p   a polynomial.
*   @param[in]  x   variable to evaluate.
*   @param[out] ppx on exit contains \f$p'(x)\f$.
*   @returns The value of \f$p(x)\f$.
*
*   @see    evalAndDerive(const Polynomial<T>& p, const U& x, U& ppx, U& pppx)
*/

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

/**
*   evaluates \f$p(x)\f$, first derivate \f$p'(x)\f$, and second derivative
*   \f$p''(x)\f$ at the same time.
*
*   Evaluates efficiently \f$p(x)\f$, first derivate \f$p'(x)\f$, and second 
*   derivative \f$p''(x)\f$ at the same time.
*
*   @param[in]  p       a polynomial.
*   @param[in]  x       variable to evaluate.
*   @param[out] ppx     on exit contains \f$p'(x)\f$.
*   @param[out] pppx    on exit contains \f$p''(x)\f$.
*   @returns            The value of \f$p(x)\f$.
*
*   @ingroup PolyEval
*   @see    evalAndDerive(const Polynomial<T>& p, const U& x, U& ppx)
*/

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

/**
*   evaluates first derivate \f$p'(x)\f$, second derivative
*   \f$p''(x)\f$, and deflates a polynomial at the same time.
*
*   Evaluates efficiently \f$p(x)\f$, first derivate \f$p'(x)\f$, and second 
*   derivative \f$p''(x)\f$, while deflating the polynomial \a p.
*
*   @param[in]  p       a polynomial.
*   @param[in]  x       variable to evaluate.
*   @param[out] ppx     on exit, contains \f$p'(x)\f$.
*   @param[out] pppx    on exit, contains \f$p''(x)\f$.
*   @param[out] q       on exit, contains \f$q(z)=\frac{p(z)}{z-x}\f$
*   @returns            The value of \f$p(x)\f$.
*
*   @ingroup PolyEval
*   @see    evalAndDerive(const Polynomial<T>& p, const U& x, U& ppx, U& ppx)
*/

template <class T, class U>
U evalDeriveAndDeflate(const Polynomial<T>& p, const U& x, U& ppx, U& pppx, Polynomial<U>& q) {
    int m = p.degree();
    if (m >= 2) {
        if (x == U(0)) {
            q.resize(m);
            std::copy(&p[1], &p[m + 1], &q[0]);
            ppx = U(p[1]);
            pppx = U(p[2] + p[2]);
            return U(p[0]);
        }
        else {
            q.resize(m);
            U px = p[m];
            ppx = U(0);
            pppx = U(0);
            for (int i = m - 1; i >= 0; --i) {
                pppx *= x; 
                pppx += ppx;
                ppx *= x; 
                ppx += px;
                q[i] = px;
                px *= x; 
                px += p[i];
            }
            pppx += pppx;
            return px;
        }
    }
    else if (m == 1) {
        q.resize(1);
        ppx = U(p[1]);
        pppx = U(0);
        q[0] = T(0);
        return (p[1] * x) + p[0];
    }
    else {
        q.resize(1);
        q[0] = U(0);
        ppx = U(0);
        pppx = U(0);
        return p[0];
    }
}

/**
*   returns optimal scale for polynomial.
*
*   Computes an optimal scale factor for polynomial \a p.  Thus ensuring that
*   evaluating the polynomial will not result in overflow or underflow.
*
*   @param[in]  p       a polynomial.
*
*   @see U scalePoly(Polynomial<T>& p, const Polynomial<U>& q);
*/

template<class T>
T getOptimalScale(const Polynomial<T>& p) {
    static FloatSpecs<T> fpSpecs;
    T hi = sqrt(fpSpecs.INFINY);
    T lo = fpSpecs.SMALNO / fpSpecs.ETA;
    T mx = 0;
    T mn = fpSpecs.INFINY;
    T x;
    for (unsigned i = 1; i < p.size(); ++i) {
        if (p[i] > mx) mx = p[i];
        if (p[i] != 0 && p[i] < mn) mn = p[i];
    }

    if (mn >= lo && mx <= hi)
        return T(1);

    x = lo / mn;
    if (x <= T(1))
        x = 1 / (sqrt(mx) * sqrt(mn));
    else
        if (fpSpecs.INFINY / x > mx)
            x = T(1);

    return pow(fpSpecs.BASE, floor(T(.5) + (log(x) / fpSpecs.LOG_BASE)));
}

/**
*   rescales a polynomial using an optimal scale.
*
*   Computes an optimal scale factor for polynomial \a p.  Thus ensuring that
*   evaluating the polynomial will not result in overflow or underflow, then
*   scales the polynomial \a p.
*
*   @param[in,out]  p       a polynomial.
*   @param[in]      q       the modulus polynomial of p.
*
*   @see    scalePoly(Polynomial<T>& p, const Polynomial<U>& q)\n
*           modulus(const Polynomial<T>& p, Polynomial<U>& q)
*/

template <class T, class U>
U scalePoly(Polynomial<T>& p, const Polynomial<U>& q) {
    U scale = getOptimalScale(q);
    p *= scale;
    return scale;
}

/**
*   computes the modulus of polynomial \a p.
*
*   Computes the modulus of polynomial \a p, such as for polynomial:
*
*   \f[p(x)=a_0+a_1{x}+a_2{x^2}+...+a_n{x^n}\f]
*
*   the resulting polynomial \a q will be :
*
*   \f[q(x)=|a_0|+|a_1|{x}+|a_2|{x^2}+...+|a_n|{x^n}\f]
*
*   @param[in]  p       a polynomial.
*   @param[out] q       the modulus polynomial p.
*/

template <class T, class U>
void modulus(const Polynomial<T>& p, Polynomial<U>& q) {
    int m = p.degree();
    q.resize(m + 1);
    for(int i = 0; i <= m; ++i) 
        q[i] = abs(p[i]);
}

/**@} Poly */
#endif //POLYNOMIAL_H