Polynomial.h

Go to the documentation of this file.

#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H

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

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

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

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

    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();
        }
    }
    void compact() {
        while((back() == T(0)) && size() > 1) pop_back();
    }
    
    int degree() const {
        int m = (int)size() - 1;
        while ((m > 0) &&( (*this)[m] == T(0))) --m;
        return m;
    }

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

    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);
        }
    }

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

    inline bool isMonic() const {
        return ((*this)[degree()] == T(1));
    }
 
    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];
    }

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

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

    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;
    }
    inline const Polynomial<T>& operator= (const std::vector<T>& q) {
        (std::vector<T>&)*this = q;
        return *this;
    }

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

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

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

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

    Polynomial<T>& operator/= (const T& s) {
        mul(*this, T(1) / s);
        return *this;
    }
    friend inline Polynomial<T> operator+ (const Polynomial<T>& p, const Polynomial<T>& q) {
        Polynomial<T> r;
        add(r, p, q);
        return r;
    }

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

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

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

    friend Polynomial<T> operator/ (const Polynomial<T>& p, const T& s) {
        Polynomial<T> r;
        mul(r, T(1) / f, p);
        return r;
    }
};

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();
}

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();
}

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();
}

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();
}

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);
    }
}

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);
    }
}

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];
    }
}

template <class T>
inline void mul(Polynomial<T>& p, const T& s) {
    mul(p, s, 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];
}

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;
    }
}

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];
        }
    }
}

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];
    }
}

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];
    }
}

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);
}

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];
    }
}

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];
    }
}

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];
    }
}

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)));
}

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

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]);
}

#endif //POLYNOMIAL_H

---