Polynomial.h

Go to the documentation of this file.
00001 #ifndef POLYNOMIAL_H
00002 #define POLYNOMIAL_H
00003 
00033 #include <vector>
00034 #include <algorithm>
00035 #include <complex>
00036 #include <stdarg.h>
00037 #include "floatspecs.h"
00038 
00085 template <class T>
00086 class Polynomial : public std::vector<T> {
00087 
00088 public:
00100     explicit Polynomial(const T& a = T(0)) {
00101         push_back(a);
00102     }
00103 
00113     explicit Polynomial(const std::vector<T>& p) : std::vector(p) { }
00114 
00127     Polynomial(unsigned n, T a0, ...) {
00128         va_list marker;
00129         int i = n;
00130         if (sizeof(T) != sizeof(float)) {
00131             push_back(a0);
00132             va_start(marker, a0);
00133             while (i-- > 0) {
00134                 T an = va_arg(marker, T);
00135                 push_back(an);
00136             }
00137             va_end(marker);
00138             compact();
00139         }
00140         else {
00141             // in c/c++, floats are pushed as doubles
00142             push_back(a0);
00143             va_start(marker, a0);
00144             while (i-- > 0) {
00145                 double an = va_arg(marker, double);
00146                 push_back(T(an));
00147             }
00148             va_end(marker);
00149             compact();
00150         }
00151     }
00168     void compact() {
00169         while((back() == T(0)) && size() > 1) pop_back();
00170     }
00171     
00188     int degree() const {
00189         int m = (int)size() - 1;
00190         while ((m > 0) &&( (*this)[m] == T(0))) --m;
00191         return m;
00192     }
00193 
00208     int degree() {
00209         compact();
00210         int m = (int)size() - 1;
00211         if (m < 0) {
00212             resize(1);
00213             m = 0;
00214         }
00215         return (int)size() - 1;
00216     }
00217 
00230     void shift(int n) {
00231         if (n < 0) {
00232             int newSize = (int)size() + n;
00233             if (newSize >= 1) {
00234                 std::copy(&(*this)[-n], &*end(), &(*this)[0]);
00235             }
00236             else {
00237                 newSize = 1;
00238                 (*this)[0] = T(0);
00239             }
00240             resize(newSize);
00241         }
00242         else if (n > 0) {
00243             resize(size() + n);
00244             std::copy(begin(), end(), begin() + n);
00245         }
00246     }
00247 
00257     void makeMonic() {
00258         compact();
00259         if (!isMonic())
00260             *this /= (*this)[degree()];
00261     }
00262 
00273     inline bool isMonic() const {
00274         return ((*this)[degree()] == T(1));
00275     }
00276  
00287     void getDerivative(Polynomial<T>& pp) const {
00288         unsigned i;
00289         T ii;
00290         pp.resize(size() - 1);
00291         for (i = 0, ii = T(1); i < pp.size(); ++i, ii += T(1))
00292             pp[i] = ii * (*this)[i + 1];
00293     }
00294 
00312     T eval(const T& x) const {
00313         return ::eval(*this, x);
00314     }
00315 
00339     T evalAndDeflate(const T& a, Polynomial<T>& q) const {
00340         return ::evalAndDeflate(*this, a, q);
00341     }
00342 
00354     T evalAndDerive(const T& x, T& ppx) const {
00355         int m = degree();
00356         T y;
00357         ppx = 0;
00358         y = (*this)[m];
00359         for (int i = m - 1; i >= 0; --i) {
00360             ppx = (ppx * x) + y;
00361             y = (y * x) + (*this)[i];
00362         }
00363         return y;
00364     }
00379     inline const Polynomial<T>& operator= (const std::vector<T>& q) {
00380         (std::vector<T>&)*this = q;
00381         return *this;
00382     }
00383 
00397     inline const Polynomial<T>& operator+= (const Polynomial<T>& q) {
00398         add(*this, q);
00399         return *this;
00400     }
00401 
00415     inline const Polynomial<T>& operator-= (const Polynomial<T>& q) {
00416         sub(*this, q);
00417         return *this;
00418     }
00419 
00433     inline const Polynomial<T> operator*= (const Polynomial<T>& p) {
00434         Polynomial<T> q = *this;
00435         mul(*this, p, q);
00436         return *this;
00437     }
00438 
00451     const Polynomial<T>& operator*= (const T& f) {
00452         mul(*this, f);
00453         return *this;
00454     }
00455 
00468     Polynomial<T>& operator/= (const T& s) {
00469         mul(*this, T(1) / s);
00470         return *this;
00471     }
00489     friend inline Polynomial<T> operator+ (const Polynomial<T>& p, const Polynomial<T>& q) {
00490         Polynomial<T> r;
00491         add(r, p, q);
00492         return r;
00493     }
00494 
00509     friend inline Polynomial<T> operator- (const Polynomial<T>& p, const Polynomial<T>& q) {
00510         Polynomial<T> r;
00511         sub(r, a, b);
00512         return r;
00513     }
00514 
00530     friend inline Polynomial<T> operator* (const Polynomial<T>& p, const Polynomial<T>& q) {
00531         Polynomial<T> r;
00532         mul(r, p, q);
00533         return r;
00534     }
00535 
00548     friend inline Polynomial<T> operator* (const T& s, const Polynomial<T>& p) {
00549         Polynomial<T> r;
00550         mul(r, s, p);
00551         return r;
00552     }
00553 
00566     friend Polynomial<T> operator/ (const Polynomial<T>& p, const T& s) {
00567         Polynomial<T> r;
00568         mul(r, T(1) / f, p);
00569         return r;
00570     }
00572 };
00573 
00589 template <class T>
00590 void add(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q) {
00591     int i, m, n;
00592     m = p.degree();
00593     n = q.degree();
00594     if (m <= n) {
00595         r.resize(n + 1);
00596         for (i = 0; i <= m; ++i)
00597             r[i] = p[i] + q[i];
00598         for ( ; i <= n; ++i)
00599             r[i] = q[i];
00600     }
00601     else {
00602         r.resize(m + 1);
00603         for (i = 0; i <= n; ++i)
00604             r[i] = p[i] + q[i];
00605         for ( ; i <= m; ++i)
00606             r[i] = p[i];
00607     }
00608     r.compact();
00609 }
00610 
00624 template <class T>
00625 void add(Polynomial<T>& p, const Polynomial<T>& q) {
00626     int m = p.degree();
00627     int n = q.degree();
00628     if (m < n)
00629         p.resize(n + 1, T(0));
00630 
00631     for (int i = 0; i <= n; ++i)
00632         p[i] += q[i];
00633     p.compact();
00634 }
00635 
00649 template <class T>
00650 void sub(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q) {
00651     int i;
00652     int m = p.degree();
00653     int n = q.degree();
00654     if (m <= n) {
00655         r.resize(n + 1);
00656         for (i = 0; i <= m; ++i)
00657             r[i] = p[i] - q[i];
00658         for ( ; i <= n; ++i)
00659             r[i] = -q[i];
00660     }
00661     else {
00662         r.resize(m + 1);
00663         for (i = 0; i <= n; ++i)
00664             r[i] = p[i] - q[i];
00665         for ( ; i <= m; ++i)
00666             r[i] = p[i];
00667     }
00668     r.compact();
00669 }
00670 
00682 template <class T>
00683 void sub(Polynomial<T>& p, const Polynomial<T>& q) {
00684     int i;
00685     int m = p.degree();
00686     int n = q.degree();
00687     if (m < n)
00688         p.resize(n + 1, T(0));
00689     for (i = 0; i <= n; ++i)
00690         p[i] -= q[i];
00691     p.compact();
00692 }
00693 
00710 template <class T>
00711 void div(Polynomial<T>& q, Polynomial<T>& r, const Polynomial<T>& u, const Polynomial<T>& v) {
00712     int na = u.degree();
00713     int nb = v.degree();
00714     if (na >= nb) {
00715         r = u;
00716         T t;
00717         q.resize(1 + na - nb);
00718         std::fill(&q[0], &q[q.size()], T(0));
00719         for (int i = na - nb; i >= 0; --i) {
00720             q[i] = r[nb + i] / v[nb];
00721             r[nb + i] = T(0);
00722             for (int j = nb + i - 1; j >= i; --j) {
00723                 r[j + 1] -= (q[i] * v[j - i]);
00724             }
00725             r[i + nb] = 0;
00726         }
00727         r.compact();
00728     }
00729     else {
00730         r = u;
00731         q.resize(1);
00732         q[0] = T(0);
00733     }
00734 }
00735 
00751 template <class T>
00752 void div(Polynomial<T>& u, Polynomial<T>& r, const Polynomial<T>& v) {
00753     int na = u.degree();
00754     int nb = v.degree();
00755     if (na >= nb) {
00756         r = u;
00757         u.resize(na - nb + 1);
00758         for (int i = na - nb; i >= 0; --i) {
00759             u[i] = r[nb + i] / v[nb];
00760             r[nb + i] = T(0);
00761             for (int j = nb + i - 1; j >= i; --j)
00762                 r[j] -= (u[i] * v[j - i]);
00763         }
00764         r.compact();
00765     }
00766     else {
00767         r = u;
00768         u.resize(1);
00769         u[0] = T(0);
00770     }
00771 }
00772 
00788 template <class T>
00789 void mul(Polynomial<T>& r, const T& s, const Polynomial<T>& p) {
00790     if (s != T(1)) {
00791         r.resize(p.degree() + 1);
00792         for (unsigned i = 0; i < r.size(); ++i)
00793             r[i] = s * p[i];
00794     }
00795 }
00796 
00810 template <class T>
00811 inline void mul(Polynomial<T>& p, const T& s) {
00812     mul(p, s, p);
00813 }
00814 
00829 template <class T>
00830 void mul(Polynomial<T>& r, const Polynomial<T>& p, const Polynomial<T>& q) {
00831     int na = p.degree();
00832     int nb = q.degree();
00833     r.resize(na + nb + 1);
00834     std::fill(&r[0], &r[r.size()], T(0));
00835     for (int i = 0; i <= na; ++i)
00836         for (int j = 0; j <= nb; ++j)
00837             r[i + j] += p[i] * q[j];
00838 }
00839 
00859 template<class T, class U>
00860 U eval(const Polynomial<T>& p, const U& x) {
00861     if (x == U(0)) 
00862         return p[0];
00863     else {
00864         int i, m = p.degree();
00865         register U px = p[m];
00866         for (i = m - 1; i >= 0; --i) {
00867             px *= x;
00868             px += p[i];
00869         }
00870         return px;
00871     }
00872 }
00873 
00895 template<class T, class U, class V>
00896 U eval(const Polynomial<T>& p, const U& x, V& e) {
00897     static FloatSpecs<V> fpSpecs;
00898     if (x == U(0)) {
00899         e = V(0);
00900         return p[0];
00901     }
00902     else {
00903         int i, m = p.degree();
00904         if (m >= 1) {
00905             U px = p[m];
00906             V mx, mpx;
00907             mx = abs(x);
00908             e = (abs(px) * fpSpecs.MRE) / (fpSpecs.ARE + fpSpecs.MRE);
00909             for (i = m - 1; i >= 0; --i) {
00910                 px *= x;
00911                 px += p[i];
00912                 //px = (px * x) + p[i];
00913                 mpx = abs(px);
00914                 e = (mx * e) + mpx;
00915             }
00916             e = (e * (fpSpecs.ARE + fpSpecs.MRE)) - (mpx * fpSpecs.MRE);
00917             return px;
00918         }
00919         else {
00920             e = V(0);
00921             return p[0];
00922         }
00923     }
00924 }
00925 
00953 template<class T, class U>
00954 U evalAndDeflate(const Polynomial<T>& p, const U& a, Polynomial<U>& q) {
00955     int i, m = p.degree();
00956     if (m > 0) {
00957         q.resize(m);
00958         if (a == T(0)) {
00959             std::copy(&p[1], &p[m + 1], &q[0]);
00960             return p[0];
00961         }
00962         else {
00963             U pa;
00964             pa = p[m];
00965             for (i = m - 1; i >= 0; --i) {
00966                 q[i] = pa;
00967                 pa *= a;
00968                 pa += p[i];
00969                 //pa = (pa * a) + p[i];
00970             }
00971             return pa;
00972         }
00973     }
00974     else {
00975         q.resize(1);
00976         q[0] = 0;
00977         return p[0];
00978     }
00979 }
00980 
01010 template<class T, class U, class V>
01011 U evalAndDeflate(const Polynomial<T>& p, const U& a, Polynomial<U>& q, V& e) {
01012     static FloatSpecs<V> fpSpecs;
01013     int i, m = p.degree();
01014     if (m > 0) {
01015         q.resize(m);
01016         if (a == T(0)) {
01017             e = 0;
01018             std::copy(&p[1], &p[m + 1], &q[0]);
01019             return p[0];
01020         }
01021         else {
01022             U pa;
01023             V mx, mpa;
01024             pa = p[m];
01025             mx = abs(a);
01026             e = (abs(pa) * fpSpecs.MRE) / (fpSpecs.ARE + fpSpecs.MRE);
01027             for (i = m - 1; i >= 0; --i) {
01028                 q[i] = pa;
01029                 pa *= a;
01030                 pa += p[i];
01031                 mpa = abs(pa);
01032                 e = (mx * e) + mpa;
01033             }
01034             e = (e * (fpSpecs.ARE + fpSpecs.MRE)) - (mpa * fpSpecs.MRE);
01035             return pa;
01036         }
01037     }
01038     else {
01039         e = 0;
01040         q.resize(1);
01041         q[0] = 0;
01042         return p[0];
01043     }
01044 }
01045 
01062 template <class T, class U>
01063 U evalError(const Polynomial<T>& p, const U& mx) {
01064     static FloatSpecs<U> fpSpecs;
01065     T px;
01066     U e, mpx;
01067     int m = p.degree();
01068     px = p[m];
01069     e = (abs(px) * fpSpecs.MRE) / (fpSpecs.ARE + fpSpecs.MRE);
01070     for (int i = m - 1; i >= 0; --i) {
01071         px *= mx;
01072         px += p[i];
01073         mpx = abs(px);
01074         e = (mx * e) + mpx;
01075     }
01076     return (e * (fpSpecs.ARE + fpSpecs.MRE)) - (mpx * fpSpecs.MRE);
01077 }
01078 
01094 template <class T, class U>
01095 U evalAndDerive(const Polynomial<T>& p, const U& x, U& ppx) {
01096     int m = p.degree();
01097     if (m >= 1) {
01098         if (x == U(0)) {
01099             ppx = U(p[1]);
01100             return U(p[0]);
01101         }
01102         else {
01103             int m = p.degree();
01104             U px = p[m];
01105             ppx = U(0);
01106             for (int i = m - 1; i >= 0; --i) {
01107                 ppx = (ppx * x) + px;
01108                 px = (px * x) + p[i];
01109                 //ppx *= x; 
01110                 //ppx += px;
01111                 //px *= x; 
01112                 //px += p[i];
01113             }
01114             return px;
01115         }
01116     }
01117     else {
01118         ppx = U(0);
01119         return p[0];
01120     }
01121 }
01122 
01140 template <class T, class U>
01141 U evalAndDerive(const Polynomial<T>& p, const U& x, U& ppx, U& pppx) {
01142     int m = p.degree();
01143     if (m >= 2) {
01144         if (x == U(0)) {
01145             ppx = U(p[1]);
01146             pppx = U(p[2] + p[2]);
01147             return U(p[0]);
01148         }
01149         else {
01150             int m = p.degree();
01151             U px = p[m];
01152             ppx = U(0);
01153             pppx = U(0);
01154             for (int i = m - 1; i >= 0; --i) {
01155                 pppx *= x; 
01156                 pppx += ppx;
01157                 ppx *= x; 
01158                 ppx += px;
01159                 px *= x; 
01160                 px += p[i];
01161             }
01162             pppx += pppx;
01163             return px;
01164         }
01165     }
01166     else if (m == 1) {
01167         ppx = U(p[1]);
01168         pppx = U(0);
01169         return (p[1] * x) + p[0];
01170     }
01171     else {
01172         ppx = U(0);
01173         pppx = U(0);
01174         return p[0];
01175     }
01176 }
01177 
01196 template <class T, class U>
01197 U evalDeriveAndDeflate(const Polynomial<T>& p, const U& x, U& ppx, U& pppx, Polynomial<U>& q) {
01198     int m = p.degree();
01199     if (m >= 2) {
01200         if (x == U(0)) {
01201             q.resize(m);
01202             std::copy(&p[1], &p[m + 1], &q[0]);
01203             ppx = U(p[1]);
01204             pppx = U(p[2] + p[2]);
01205             return U(p[0]);
01206         }
01207         else {
01208             q.resize(m);
01209             U px = p[m];
01210             ppx = U(0);
01211             pppx = U(0);
01212             for (int i = m - 1; i >= 0; --i) {
01213                 pppx *= x; 
01214                 pppx += ppx;
01215                 ppx *= x; 
01216                 ppx += px;
01217                 q[i] = px;
01218                 px *= x; 
01219                 px += p[i];
01220             }
01221             pppx += pppx;
01222             return px;
01223         }
01224     }
01225     else if (m == 1) {
01226         q.resize(1);
01227         ppx = U(p[1]);
01228         pppx = U(0);
01229         q[0] = T(0);
01230         return (p[1] * x) + p[0];
01231     }
01232     else {
01233         q.resize(1);
01234         q[0] = U(0);
01235         ppx = U(0);
01236         pppx = U(0);
01237         return p[0];
01238     }
01239 }
01240 
01252 template<class T>
01253 T getOptimalScale(const Polynomial<T>& p) {
01254     static FloatSpecs<T> fpSpecs;
01255     T hi = sqrt(fpSpecs.INFINY);
01256     T lo = fpSpecs.SMALNO / fpSpecs.ETA;
01257     T mx = 0;
01258     T mn = fpSpecs.INFINY;
01259     T x;
01260     for (unsigned i = 1; i < p.size(); ++i) {
01261         if (p[i] > mx) mx = p[i];
01262         if (p[i] != 0 && p[i] < mn) mn = p[i];
01263     }
01264 
01265     if (mn >= lo && mx <= hi)
01266         return T(1);
01267 
01268     x = lo / mn;
01269     if (x <= T(1))
01270         x = 1 / (sqrt(mx) * sqrt(mn));
01271     else
01272         if (fpSpecs.INFINY / x > mx)
01273             x = T(1);
01274 
01275     return pow(fpSpecs.BASE, floor(T(.5) + (log(x) / fpSpecs.LOG_BASE)));
01276 }
01277 
01292 template <class T, class U>
01293 U scalePoly(Polynomial<T>& p, const Polynomial<U>& q) {
01294     U scale = getOptimalScale(q);
01295     p *= scale;
01296     return scale;
01297 }
01298 
01314 template <class T, class U>
01315 void modulus(const Polynomial<T>& p, Polynomial<U>& q) {
01316     int m = p.degree();
01317     q.resize(m + 1);
01318     for(int i = 0; i <= m; ++i) 
01319         q[i] = abs(p[i]);
01320 }
01321 
01323 #endif //POLYNOMIAL_H

Generated on Mon Aug 21 21:00:38 2006 for The Polynomials Templates Library by  doxygen 1.4.5