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
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
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
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
01110
01111
01112
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