Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_Tay_TaylorImp.hpp
Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 #include "Sacado_ConfigDefs.h"
00033 #include "Sacado_DynamicArrayTraits.hpp"
00034 #include <ostream>  // for std::ostream
00035 
00036 namespace Sacado {
00037 namespace Tay {
00038 
00039 template <typename T> 
00040 Taylor<T>::TaylorData::
00041 TaylorData() :
00042   coeff_(NULL), deg_(-1), len_(0) 
00043 {
00044 }
00045 
00046 template <typename T> 
00047 Taylor<T>::TaylorData::
00048 TaylorData(const T& x) :
00049   coeff_(), deg_(0), len_(1) 
00050 {
00051   coeff_ = Sacado::ds_array<T>::get_and_fill(len_);
00052   coeff_[0] = x;
00053 }
00054 
00055 template <typename T> 
00056 Taylor<T>::TaylorData::
00057 TaylorData(unsigned int d, const T& x) :
00058   coeff_(), deg_(d), len_(d+1) 
00059 {
00060   coeff_ = Sacado::ds_array<T>::get_and_fill(len_);
00061   coeff_[0] = x;
00062 }
00063 
00064 template <typename T> 
00065 Taylor<T>::TaylorData::
00066 TaylorData(unsigned int d) :
00067   coeff_(), deg_(d), len_(d+1) 
00068 {
00069   coeff_ = Sacado::ds_array<T>::get_and_fill(len_);
00070 }
00071 
00072 template <typename T> 
00073 Taylor<T>::TaylorData::
00074 TaylorData(unsigned int d, unsigned int l) :
00075   coeff_(), deg_(d), len_(l) 
00076 {
00077   coeff_ = Sacado::ds_array<T>::get_and_fill(len_);
00078 }
00079 
00080 template <typename T> 
00081 Taylor<T>::TaylorData::
00082 TaylorData(const typename Taylor<T>::TaylorData& x) :
00083   coeff_(), deg_(x.deg_), len_(x.deg_+1) 
00084 {
00085   coeff_ = Sacado::ds_array<T>::get_and_fill(x.coeff_, len_);
00086 }
00087 
00088 template <typename T> 
00089 Taylor<T>::TaylorData::
00090 ~TaylorData()
00091 {
00092   if (len_ > 0)
00093     Sacado::ds_array<T>::destroy_and_release(coeff_, len_);
00094 }
00095 
00096 template <typename T> 
00097 typename Taylor<T>::TaylorData&
00098 Taylor<T>::TaylorData::
00099 operator=(const typename Taylor<T>::TaylorData& x) 
00100 {
00101   if (len_ < x.deg_+1) {
00102     Sacado::ds_array<T>::destroy_and_release(coeff_, len_);
00103     len_ = x.deg_+1;
00104     deg_ = x.deg_;
00105     coeff_ = Sacado::ds_array<T>::get_and_fill(x.coeff_, len_);
00106   }
00107   else {
00108     deg_ = x.deg_;
00109     Sacado::ds_array<T>::copy(x.coeff_, coeff_, deg_+1);
00110   }
00111   
00112   return *this;
00113 }
00114 
00115 template <typename T> 
00116 Taylor<T>::
00117 Taylor() :
00118   th(new TaylorData)
00119 {
00120 }
00121 
00122 template <typename T> 
00123 Taylor<T>::
00124 Taylor(const T& x) :
00125   th(new TaylorData(x))
00126 {
00127 }
00128 
00129 template <typename T> 
00130 Taylor<T>::
00131 Taylor(const typename dummy<value_type,scalar_type>::type& x) :
00132   th(new TaylorData(value_type(x)))
00133 {
00134 }
00135 
00136 template <typename T> 
00137 Taylor<T>::
00138 Taylor(unsigned int d, const T& x) :
00139   th(new TaylorData(d, x))
00140 {
00141 }
00142 
00143 template <typename T> 
00144 Taylor<T>::
00145 Taylor(unsigned int d) :
00146   th(new TaylorData(d))
00147 {
00148 }
00149 
00150 template <typename T> 
00151 Taylor<T>::
00152 Taylor(const Taylor<T>& x) :
00153   th(x.th)
00154 {
00155 }
00156 
00157 template <typename T> 
00158 Taylor<T>::
00159 ~Taylor()
00160 {
00161 }
00162 
00163 template <typename T> 
00164 void
00165 Taylor<T>::
00166 resize(unsigned int d, bool keep_coeffs)
00167 {
00168   if (d+1 > length()) {
00169     Sacado::Handle<TaylorData> h(new TaylorData(d));
00170     if (keep_coeffs)
00171       Sacado::ds_array<T>::copy(th->coeff_, h->coeff_, th->deg_+1);
00172     th = h;
00173   }
00174   else {
00175     th.makeOwnCopy();
00176     if (!keep_coeffs)
00177       Sacado::ds_array<T>::zero(coeff(), degree()+1);
00178     th->deg_ = d;
00179   }
00180 }
00181 
00182 template <typename T> 
00183 void
00184 Taylor<T>::
00185 reserve(unsigned int d)
00186 {
00187   if (d+1 > length()) {
00188     Sacado::Handle<TaylorData> h(new TaylorData(th->deg_,d+1));
00189     Sacado::ds_array<T>::copy(th->coeff_, h->coeff_, th->deg_+1);
00190     th = h;
00191   }
00192 }
00193 
00194 template <typename T> 
00195 Taylor<T>& 
00196 Taylor<T>::
00197 operator=(const T& v) 
00198 {
00199   th.makeOwnCopy();
00200 
00201   if (th->len_ == 0) {
00202     th->len_ = 1;
00203     th->deg_ = 0;
00204     th->coeff_ = Sacado::ds_array<T>::get_and_fill(th->len_);
00205   }
00206 
00207   th->coeff_[0] = v;
00208   Sacado::ds_array<T>::zero(th->coeff_+1, th->deg_);
00209 
00210   return *this;
00211 }
00212 
00213 template <typename T> 
00214 Taylor<T>& 
00215 Taylor<T>::
00216 operator=(const typename dummy<value_type,scalar_type>::type& v) 
00217 {
00218   return operator=(value_type(v));
00219 }
00220 
00221 template <typename T> 
00222 Taylor<T>& 
00223 Taylor<T>::
00224 operator=(const Taylor<T>& x) 
00225 {
00226   th = x.th;
00227   return *this;
00228 }
00229 
00230 template <typename T> 
00231 Taylor<T>
00232 Taylor<T>::
00233 operator+() const
00234 {
00235   return *this;
00236 }
00237 
00238 template <typename T> 
00239 Taylor<T> 
00240 Taylor<T>::
00241 operator-() const
00242 {
00243   Taylor<T> x;
00244   x.th->deg_ = th->deg_;
00245   x.th->len_ = th->deg_+1;
00246   x.th->coeff_ = Sacado::ds_array<T>::get_and_fill(x.th->len_);
00247   for (unsigned int i=0; i<=th->deg_; i++)
00248     x.th->coeff_[i] = -th->coeff_[i];
00249 
00250   return x;
00251 }
00252 
00253 template <typename T> 
00254  Taylor<T>& 
00255 Taylor<T>::
00256 operator+=(const T& v)
00257 {
00258   th.makeOwnCopy();
00259 
00260   th->coeff_[0] += v;
00261 
00262   return *this;
00263 }
00264 
00265 template <typename T> 
00266 Taylor<T>& 
00267 Taylor<T>::
00268 operator-=(const T& v)
00269 {
00270   th.makeOwnCopy();
00271 
00272   th->coeff_[0] -= v;
00273 
00274   return *this;
00275 }
00276 
00277 template <typename T> 
00278 Taylor<T>& 
00279 Taylor<T>::
00280 operator*=(const T& v)
00281 {
00282   th.makeOwnCopy();
00283 
00284   for (unsigned int i=0; i<=th->deg_; i++)
00285     th->coeff_[i] *= v;
00286 
00287   return *this;
00288 }
00289 
00290 template <typename T> 
00291 Taylor<T>& 
00292 Taylor<T>::
00293 operator/=(const T& v)
00294 {
00295   th.makeOwnCopy();
00296 
00297   for (unsigned int i=0; i<=th->deg_; i++)
00298     th->coeff_[i] /= v;
00299 
00300   return *this;
00301 }
00302 
00303 template <typename T> 
00304 Taylor<T>& 
00305 Taylor<T>::
00306 operator+=(const Taylor<T>& x)
00307 {
00308   th.makeOwnCopy();
00309 
00310   unsigned int d = degree();
00311   unsigned int xd = x.degree();
00312   unsigned int dmin = xd < d ? xd : d;
00313 
00314   unsigned int l = xd + 1;
00315   bool need_resize = l > length();
00316   T* c;
00317   if (need_resize) {
00318     c = Sacado::ds_array<T>::get_and_fill(l);
00319     for (unsigned int i=0; i<=d; i++)
00320       c[i] = fastAccessCoeff(i);
00321   }
00322   else
00323     c = th->coeff_;
00324   T* xc = x.th->coeff_;
00325 
00326   for (unsigned int i=0; i<=dmin; i++)
00327     c[i] += xc[i];
00328   if (th->deg_ < xd) {
00329     for (unsigned int i=d+1; i<=xd; i++)
00330       c[i] = xc[i];
00331     th->deg_ = xd;
00332   }
00333 
00334   if (need_resize) {
00335     Sacado::ds_array<T>::destroy_and_release(th->coeff_, th->len_);
00336     th->len_ = l;
00337     th->coeff_ = c;
00338   }
00339 
00340   return *this;
00341 }
00342 
00343 template <typename T> 
00344 Taylor<T>& 
00345 Taylor<T>::
00346 operator-=(const Taylor<T>& x)
00347 {
00348   th.makeOwnCopy();
00349 
00350   unsigned int d = degree();
00351   unsigned int xd = x.degree();
00352   unsigned int dmin = xd < d ? xd : d;
00353 
00354   unsigned int l = xd + 1;
00355   bool need_resize = l > length();
00356   T* c;
00357   if (need_resize) {
00358     c = Sacado::ds_array<T>::get_and_fill(l);
00359     for (unsigned int i=0; i<=d; i++)
00360       c[i] = fastAccessCoeff(i);
00361   }
00362   else
00363     c = th->coeff_;
00364   T* xc = x.th->coeff_;
00365 
00366   for (unsigned int i=0; i<=dmin; i++)
00367     c[i] -= xc[i];
00368   if (d < xd) {
00369     for (unsigned int i=d+1; i<=xd; i++)
00370       c[i] = -xc[i];
00371     th->deg_ = xd;
00372   }
00373 
00374   if (need_resize) {
00375     Sacado::ds_array<T>::destroy_and_release(th->coeff_, th->len_);
00376     th->len_ = l;
00377     th->coeff_ = c;
00378   }
00379 
00380   return *this;
00381 }
00382 
00383 template <typename T> 
00384 Taylor<T>& 
00385 Taylor<T>::
00386 operator*=(const Taylor<T>& x)
00387 {
00388   unsigned int d = degree();
00389   unsigned int xd = x.degree();
00390 
00391 #ifdef SACADO_DEBUG
00392   if ((xd != d) && (xd != 0) && (d != 0))
00393     throw "Taylor Error:  Attempt to assign with incompatible degrees";
00394 #endif
00395   
00396   T* c = th->coeff_; 
00397   T* xc = x.th->coeff_;
00398   
00399   if (d > 0 && xd > 0) {
00400     Sacado::Handle<TaylorData> h(new TaylorData(d));
00401     T* cc = h->coeff_;
00402     T tmp;
00403     for (unsigned int i=0; i<=d; i++) {
00404       tmp = T(0.0);
00405       for (unsigned int k=0; k<=i; ++k)
00406   tmp += c[k]*xc[i-k];
00407       cc[i] = tmp;
00408     }
00409     th = h;
00410   }
00411   else if (d > 0) {
00412     th.makeOwnCopy();
00413     for (unsigned int i=0; i<=d; i++)
00414       c[i] = c[i]*xc[0];
00415   }
00416   else if (xd >= 0) {
00417     if (length() < xd+1) {
00418       Sacado::Handle<TaylorData> h(new TaylorData(xd));
00419       T* cc = h->coeff_;
00420       for (unsigned int i=0; i<=xd; i++)
00421   cc[i] = c[0]*xc[i];
00422       th = h;
00423     }
00424     else {
00425       th.makeOwnCopy();
00426       for (unsigned int i=d; i>=0; i--)
00427   c[i] = c[0]*xc[i];
00428     }
00429   }
00430 
00431   return *this;
00432 }
00433 
00434 template <typename T> 
00435 Taylor<T>& 
00436 Taylor<T>::
00437 operator/=(const Taylor<T>& x)
00438 {
00439   unsigned int d = degree();
00440   unsigned int xd = x.degree();
00441 
00442 #ifdef SACADO_DEBUG
00443   if ((xd != d) && (xd != 0) && (d != 0))
00444     throw "Taylor Error:  Attempt to assign with incompatible degrees";
00445 #endif
00446   
00447   T* c = th->coeff_; 
00448   T* xc = x.th->coeff_;
00449   
00450   if (d > 0 && xd > 0) {
00451     Sacado::Handle<TaylorData> h(new TaylorData(d));
00452     T* cc = h->coeff_;
00453     T tmp;
00454     for(unsigned int i=0; i<=d; i++) {
00455       tmp = c[i];
00456       for (unsigned int k=1; k<=i; ++k)
00457   tmp -= xc[k]*cc[i-k];
00458       cc[i] = tmp / xc[0];
00459     }
00460     th = h;
00461   }
00462   else if (d > 0) {
00463     th.makeOwnCopy();
00464     for (unsigned int i=0; i<=d; i++)
00465       c[i] = c[i]/xc[0];
00466   }
00467   else if (xd >= 0) {
00468     Sacado::Handle<TaylorData> h(new TaylorData(xd));
00469     T* cc = h->coeff_;
00470     T tmp;
00471     cc[0] = c[0] / xc[0];
00472     for (unsigned int i=1; i<=xd; i++) {
00473       tmp = T(0.0);
00474       for (unsigned int k=1; k<=i; ++k)
00475   tmp -= xc[k]*cc[i-k];
00476       cc[i] = tmp / xc[0];
00477     }
00478     th = h;
00479   }
00480 
00481   return *this;
00482 }
00483 
00484 template <typename T>
00485 void
00486 Taylor<T>::
00487 resizeCoeffs(unsigned int len)
00488 {
00489   if (th->coeff_)
00490     Sacado::ds_array<T>::destroy_and_release(th->coeff_, th->len_);
00491   th->len_ = len;
00492   th->coeff_ = Sacado::ds_array<T>::get_and_fill(th->len_);
00493 }
00494 
00495 template <typename T>
00496 Taylor<T>
00497 operator+(const Taylor<T>& a, 
00498     const Taylor<T>& b)
00499 {
00500   unsigned int da = a.degree();
00501   unsigned int db = b.degree();
00502   unsigned int dc = da > db ? da : db;
00503 
00504 #ifdef SACADO_DEBUG
00505   if ((da != db) && (da != 0) && (db != 0))
00506     throw "operator+():  Arguments have incompatible degrees!";
00507 #endif
00508 
00509   Taylor<T> c(dc);
00510   const T* ca = a.coeff();
00511   const T* cb = b.coeff();
00512   T* cc = c.coeff();
00513 
00514   if (da > 0 && db > 0) {
00515     for (unsigned int i=0; i<=dc; i++)
00516       cc[i] = ca[i] + cb[i];
00517   }
00518   else if (da > 0) {
00519     cc[0] = ca[0] + cb[0];
00520     for (unsigned int i=1; i<=dc; i++)
00521       cc[i] = ca[i];
00522   }
00523   else if (db >= 0) {
00524     cc[0] = ca[0] + cb[0];
00525     for (unsigned int i=1; i<=dc; i++)
00526       cc[i] = cb[i];
00527   }
00528 
00529   return c;
00530 }
00531 
00532 template <typename T>
00533 Taylor<T>
00534 operator+(const T& a, const Taylor<T>& b)
00535 {
00536   unsigned int dc = b.degree();
00537 
00538   Taylor<T> c(dc);
00539   const T* cb = b.coeff();
00540   T* cc = c.coeff();
00541 
00542   cc[0] = a + cb[0];
00543   for (unsigned int i=1; i<=dc; i++)
00544     cc[i] = cb[i];
00545 
00546   return c;
00547 }
00548 
00549 template <typename T>
00550 Taylor<T>
00551 operator+(const Taylor<T>& a, const T& b)
00552 {
00553   unsigned int dc = a.degree();
00554 
00555   Taylor<T> c(dc);
00556   const T* ca = a.coeff();
00557   T* cc = c.coeff();
00558 
00559   cc[0] = ca[0] + b;
00560   for (unsigned int i=1; i<=dc; i++)
00561     cc[i] = ca[i];
00562 
00563   return c;
00564 }
00565 
00566 template <typename T>
00567 Taylor<T>
00568 operator-(const Taylor<T>& a, 
00569     const Taylor<T>& b)
00570 {
00571   unsigned int da = a.degree();
00572   unsigned int db = b.degree();
00573   unsigned int dc = da > db ? da : db;
00574 
00575 #ifdef SACADO_DEBUG
00576   if ((da != db) && (da != 0) && (db != 0))
00577     throw "operator+():  Arguments have incompatible degrees!";
00578 #endif
00579 
00580   Taylor<T> c(dc);
00581   const T* ca = a.coeff();
00582   const T* cb = b.coeff();
00583   T* cc = c.coeff();
00584 
00585   if (da > 0 && db > 0) {
00586     for (unsigned int i=0; i<=dc; i++)
00587       cc[i] = ca[i] - cb[i];
00588   }
00589   else if (da > 0) {
00590     cc[0] = ca[0] - cb[0];
00591     for (unsigned int i=1; i<=dc; i++)
00592       cc[i] = ca[i];
00593   }
00594   else if (db >= 0) {
00595     cc[0] = ca[0] - cb[0];
00596     for (unsigned int i=1; i<=dc; i++)
00597       cc[i] = -cb[i];
00598   }
00599 
00600   return c;
00601 }
00602 
00603 template <typename T>
00604 Taylor<T>
00605 operator-(const T& a, const Taylor<T>& b)
00606 {
00607   unsigned int dc = b.degree();
00608 
00609   Taylor<T> c(dc);
00610   const T* cb = b.coeff();
00611   T* cc = c.coeff();
00612 
00613   cc[0] = a - cb[0];
00614   for (unsigned int i=1; i<=dc; i++)
00615     cc[i] = -cb[i];
00616 
00617   return c;
00618 }
00619 
00620 template <typename T>
00621 Taylor<T>
00622 operator-(const Taylor<T>& a, const T& b)
00623 {
00624   unsigned int dc = a.degree();
00625 
00626   Taylor<T> c(dc);
00627   const T* ca = a.coeff();
00628   T* cc = c.coeff();
00629 
00630   cc[0] = ca[0] - b;
00631   for (unsigned int i=1; i<=dc; i++)
00632     cc[i] = ca[i];
00633 
00634   return c;
00635 }
00636 
00637 template <typename T>
00638 Taylor<T>
00639 operator*(const Taylor<T>& a, 
00640     const Taylor<T>& b)
00641 {
00642   unsigned int da = a.degree();
00643   unsigned int db = b.degree();
00644   unsigned int dc = da > db ? da : db;
00645 
00646 #ifdef SACADO_DEBUG
00647   if ((da != db) && (da != 0) && (db != 0))
00648     throw "operator+():  Arguments have incompatible degrees!";
00649 #endif
00650 
00651   Taylor<T> c(dc);
00652   const T* ca = a.coeff();
00653   const T* cb = b.coeff();
00654   T* cc = c.coeff();
00655 
00656   if (da > 0 && db > 0) {
00657     T tmp;
00658     for (unsigned int i=0; i<=dc; i++) {
00659       tmp = T(0.0);
00660       for (unsigned int k=0; k<=i; k++)
00661   tmp += ca[k]*cb[i-k];
00662       cc[i] = tmp;
00663     }
00664   }
00665   else if (da > 0) {
00666     for (unsigned int i=0; i<=dc; i++)
00667       cc[i] = ca[i]*cb[0];
00668   }
00669   else if (db >= 0) {
00670     for (unsigned int i=0; i<=dc; i++)
00671       cc[i] = ca[0]*cb[i];
00672   }
00673 
00674   return c;
00675 }
00676 
00677 template <typename T>
00678 Taylor<T>
00679 operator*(const T& a, const Taylor<T>& b)
00680 {
00681   unsigned int dc = b.degree();
00682 
00683   Taylor<T> c(dc);
00684   const T* cb = b.coeff();
00685   T* cc = c.coeff();
00686 
00687   for (unsigned int i=0; i<=dc; i++)
00688     cc[i] = a*cb[i];
00689 
00690   return c;
00691 }
00692 
00693 template <typename T>
00694 Taylor<T>
00695 operator*(const Taylor<T>& a, const T& b)
00696 {
00697   unsigned int dc = a.degree();
00698 
00699   Taylor<T> c(dc);
00700   const T* ca = a.coeff();
00701   T* cc = c.coeff();
00702 
00703   for (unsigned int i=0; i<=dc; i++)
00704     cc[i] = ca[i]*b;
00705 
00706   return c;
00707 }
00708 
00709 template <typename T>
00710 Taylor<T>
00711 operator/(const Taylor<T>& a, 
00712     const Taylor<T>& b)
00713 {
00714   unsigned int da = a.degree();
00715   unsigned int db = b.degree();
00716   unsigned int dc = da > db ? da : db;
00717 
00718 #ifdef SACADO_DEBUG
00719   if ((da != db) && (da != 0) && (db != 0))
00720     throw "operator+():  Arguments have incompatible degrees!";
00721 #endif
00722 
00723   Taylor<T> c(dc);
00724   const T* ca = a.coeff();
00725   const T* cb = b.coeff();
00726   T* cc = c.coeff();
00727 
00728   if (da > 0 && db > 0) {
00729     T tmp;
00730     for (unsigned int i=0; i<=dc; i++) {
00731       tmp = ca[i];
00732       for (unsigned int k=0; k<=i; k++)
00733   tmp -= cb[k]*cc[i-k];
00734       cc[i] = tmp / cb[0];
00735     }
00736   }
00737   else if (da > 0) {
00738     for (unsigned int i=0; i<=dc; i++)
00739       cc[i] = ca[i]/cb[0];
00740   }
00741   else if (db >= 0) {
00742     T tmp;
00743     cc[0] = ca[0] / cb[0];
00744     for (unsigned int i=1; i<=dc; i++) {
00745       tmp = T(0.0);
00746       for (unsigned int k=0; k<=i; k++)
00747   tmp -= cb[k]*cc[i-k];
00748       cc[i] = tmp / cb[0];
00749     }
00750   }
00751 
00752   return c;
00753 }
00754 
00755 template <typename T>
00756 Taylor<T>
00757 operator/(const T& a, const Taylor<T>& b)
00758 {
00759   unsigned int dc = b.degree();
00760 
00761   Taylor<T> c(dc);
00762   const T* cb = b.coeff();
00763   T* cc = c.coeff();
00764 
00765   T tmp;
00766   cc[0] = a / cb[0];
00767   for (unsigned int i=1; i<=dc; i++) {
00768     tmp = T(0.0);
00769     for (unsigned int k=0; k<=i; k++)
00770       tmp -= cb[k]*cc[i-k];
00771     cc[i] = tmp / cb[0];
00772   }
00773 
00774   return c;
00775 }
00776 
00777 template <typename T>
00778 Taylor<T>
00779 operator/(const Taylor<T>& a, const T& b)
00780 {
00781   unsigned int dc = a.degree();
00782 
00783   Taylor<T> c(dc);
00784   const T* ca = a.coeff();
00785   T* cc = c.coeff();
00786 
00787   for (unsigned int i=0; i<=dc; i++)
00788     cc[i] = ca[i]/b;
00789 
00790   return c;
00791 }
00792 
00793 template <typename T>
00794 Taylor<T>
00795 exp(const Taylor<T>& a)
00796 {
00797   unsigned int dc = a.degree();
00798 
00799   Taylor<T> c(dc);
00800   const T* ca = a.coeff();
00801   T* cc = c.coeff();
00802 
00803   T tmp;
00804   cc[0] = std::exp(ca[0]);
00805   for (unsigned int i=1; i<=dc; i++) {
00806     tmp = T(0.0);
00807     for (unsigned int k=1; k<=i; k++)
00808       tmp += k*cc[i-k]*ca[k];
00809     cc[i] = tmp / i;
00810   }
00811 
00812   return c;
00813 }
00814 
00815 template <typename T>
00816 Taylor<T>
00817 log(const Taylor<T>& a)
00818 {
00819   unsigned int dc = a.degree();
00820 
00821   Taylor<T> c(dc);
00822   const T* ca = a.coeff();
00823   T* cc = c.coeff();
00824 
00825   T tmp;
00826   cc[0] = std::log(ca[0]);
00827   for (unsigned int i=1; i<=dc; i++) {
00828     tmp = i*ca[i];
00829     for (unsigned int k=1; k<=i-1; k++)
00830       tmp -= k*ca[i-k]*cc[k];
00831     cc[i] = tmp / (i*ca[0]);
00832   }
00833 
00834   return c;
00835 }
00836 
00837 template <typename T>
00838 Taylor<T>
00839 log10(const Taylor<T>& a)
00840 {
00841   return log(a) / std::log(10.0);
00842 }
00843 
00844 template <typename T>
00845 Taylor<T>
00846 sqrt(const Taylor<T>& a)
00847 {
00848   unsigned int dc = a.degree();
00849 
00850   Taylor<T> c(dc);
00851   const T* ca = a.coeff();
00852   T* cc = c.coeff();
00853 
00854   T tmp;
00855   cc[0] = std::sqrt(ca[0]);
00856   for (unsigned int i=1; i<=dc; i++) {
00857     tmp = ca[i];
00858     for (unsigned int k=1; k<=i-1; k++)
00859       tmp -= cc[k]*cc[i-k];
00860     cc[i] = tmp / (2.0*cc[0]);
00861   }
00862 
00863   return c;
00864 }
00865 
00866 template <typename T>
00867 Taylor<T>
00868 pow(const Taylor<T>& a,
00869     const Taylor<T>& b)
00870 {
00871   return exp(b*log(a));
00872 }
00873 
00874 template <typename T>
00875 Taylor<T>
00876 pow(const T& a,
00877     const Taylor<T>& b)
00878 {
00879   return exp(b*std::log(a));
00880 }
00881 
00882 template <typename T>
00883 Taylor<T>
00884 pow(const Taylor<T>& a,
00885     const T& b)
00886 {
00887   return exp(b*log(a));
00888 }
00889 
00890 template <typename T>
00891 void
00892 sincos(const Taylor<T>& a,
00893        Taylor<T>& s,
00894        Taylor<T>& c)
00895 {
00896   unsigned int dc = a.degree();
00897   if (s.degree() != dc)
00898     s.resize(dc, false);
00899   if (c.degree() != dc)
00900     c.resize(dc, false);
00901 
00902   const T* ca = a.coeff();
00903   T* cs = s.coeff();
00904   T* cc = c.coeff();
00905 
00906   T tmp1;
00907   T tmp2;
00908   cs[0] = std::sin(ca[0]);
00909   cc[0] = std::cos(ca[0]);
00910   for (unsigned int i=1; i<=dc; i++) {
00911     tmp1 = T(0.0);
00912     tmp2 = T(0.0);
00913     for (unsigned int k=1; k<=i; k++) {
00914       tmp1 += k*ca[k]*cc[i-k];
00915       tmp2 -= k*ca[k]*cs[i-k];
00916     }
00917     cs[i] = tmp1 / i;
00918     cc[i] = tmp2 / i;
00919   }
00920 }
00921 
00922 template <typename T>
00923 Taylor<T>
00924 sin(const Taylor<T>& a)
00925 {
00926   unsigned int dc = a.degree();
00927   Taylor<T> s(dc);
00928   Taylor<T> c(dc);
00929   sincos(a, s, c);
00930 
00931   return s;
00932 }
00933 
00934 template <typename T>
00935 Taylor<T>
00936 cos(const Taylor<T>& a)
00937 {
00938   unsigned int dc = a.degree();
00939   Taylor<T> s(dc);
00940   Taylor<T> c(dc);
00941   sincos(a, s, c);
00942 
00943   return c;
00944 }
00945 
00946 template <typename T>
00947 Taylor<T>
00948 tan(const Taylor<T>& a)
00949 {
00950   unsigned int dc = a.degree();
00951   Taylor<T> s(dc);
00952   Taylor<T> c(dc);
00953   
00954   sincos(a, s, c);
00955 
00956   return s / c;
00957 }
00958 
00959 template <typename T>
00960 void
00961 sinhcosh(const Taylor<T>& a,
00962    Taylor<T>& s,
00963    Taylor<T>& c)
00964 {
00965   unsigned int dc = a.degree();
00966   if (s.degree() != dc)
00967     s.resize(dc, false);
00968   if (c.degree() != dc)
00969     c.resize(dc, false);
00970 
00971   const T* ca = a.coeff();
00972   T* cs = s.coeff();
00973   T* cc = c.coeff();
00974 
00975   T tmp1;
00976   T tmp2;
00977   cs[0] = std::sinh(ca[0]);
00978   cc[0] = std::cosh(ca[0]);
00979   for (unsigned int i=1; i<=dc; i++) {
00980     tmp1 = T(0.0);
00981     tmp2 = T(0.0);
00982     for (unsigned int k=1; k<=i; k++) {
00983       tmp1 += k*ca[k]*cc[i-k];
00984       tmp2 += k*ca[k]*cs[i-k];
00985     }
00986     cs[i] = tmp1 / i;
00987     cc[i] = tmp2 / i;
00988   }
00989 }
00990 
00991 template <typename T>
00992 Taylor<T>
00993 sinh(const Taylor<T>& a)
00994 {
00995   unsigned int dc = a.degree();
00996   Taylor<T> s(dc);
00997   Taylor<T> c(dc);
00998   sinhcosh(a, s, c);
00999 
01000   return s;
01001 }
01002 
01003 template <typename T>
01004 Taylor<T>
01005 cosh(const Taylor<T>& a)
01006 {
01007   unsigned int dc = a.degree();
01008   Taylor<T> s(dc);
01009   Taylor<T> c(dc);
01010   sinhcosh(a, s, c);
01011 
01012   return c;
01013 }
01014 
01015 template <typename T>
01016 Taylor<T>
01017 tanh(const Taylor<T>& a)
01018 {
01019   unsigned int dc = a.degree();
01020   Taylor<T> s(dc);
01021   Taylor<T> c(dc);
01022   
01023   sinhcosh(a, s, c);
01024 
01025   return s / c;
01026 }
01027 
01028 template <typename T>
01029 Taylor<T>
01030 quad(const T& c0,
01031      const Taylor<T>& a,
01032      const Taylor<T>& b)
01033 {
01034   unsigned int dc = a.degree();
01035 
01036   Taylor<T> c(dc);
01037   const T* ca = a.coeff();
01038   const T* cb = b.coeff();
01039   T* cc = c.coeff();
01040 
01041   T tmp;
01042   cc[0] = c0;
01043   for (unsigned int i=1; i<=dc; i++) {
01044     tmp = T(0.0);
01045     for (unsigned int k=1; k<=i; k++) 
01046       tmp += k*ca[k]*cb[i-k];
01047     cc[i] = tmp / i;
01048   }
01049 
01050   return c;
01051 }
01052 
01053 template <typename T>
01054 Taylor<T>
01055 acos(const Taylor<T>& a)
01056 {
01057   Taylor<T> b = -1.0 / sqrt(1.0 - a*a);
01058   return quad(std::acos(a.coeff(0)), a, b);
01059 }
01060 
01061 template <typename T>
01062 Taylor<T>
01063 asin(const Taylor<T>& a)
01064 {
01065   Taylor<T> b = 1.0 / sqrt(1.0 - a*a);
01066   return quad(std::asin(a.coeff(0)), a, b);
01067 }
01068 
01069 template <typename T>
01070 Taylor<T>
01071 atan(const Taylor<T>& a)
01072 {
01073   Taylor<T> b = 1.0 / (1.0 + a*a);
01074   return quad(std::atan(a.coeff(0)), a, b);
01075 }
01076 
01077 template <typename T>
01078 Taylor<T>
01079 atan2(const Taylor<T>& a,
01080       const Taylor<T>& b)
01081 {
01082   Taylor<T> c = atan(a/b);
01083   c.fastAccessCoeff(0) = atan2(a.coeff(0),b.coeff(0));
01084 }
01085 
01086 template <typename T>
01087 Taylor<T>
01088 atan2(const T& a,
01089       const Taylor<T>& b)
01090 {
01091   Taylor<T> c = atan(a/b);
01092   c.fastAccessCoeff(0) = atan2(a,b.coeff(0));
01093 }
01094 
01095 template <typename T>
01096 Taylor<T>
01097 atan2(const Taylor<T>& a,
01098       const T& b)
01099 {
01100   Taylor<T> c = atan(a/b);
01101   c.fastAccessCoeff(0) = atan2(a.coeff(0),b);
01102 }
01103 
01104 template <typename T>
01105 Taylor<T>
01106 acosh(const Taylor<T>& a)
01107 {
01108   Taylor<T> b = -1.0 / sqrt(1.0 - a*a);
01109   return quad(acosh(a.coeff(0)), a, b);
01110 }
01111 
01112 template <typename T>
01113 Taylor<T>
01114 asinh(const Taylor<T>& a)
01115 {
01116   Taylor<T> b = 1.0 / sqrt(a*a - 1.0);
01117   return quad(asinh(a.coeff(0)), a, b);
01118 }
01119 
01120 template <typename T>
01121 Taylor<T>
01122 atanh(const Taylor<T>& a)
01123 {
01124   Taylor<T> b = 1.0 / (1.0 - a*a);
01125   return quad(atanh(a.coeff(0)), a, b);
01126 }
01127 
01128 template <typename T>
01129 Taylor<T>
01130 fabs(const Taylor<T>& a)
01131 {
01132   if (a.coeff(0) >= 0)
01133     return a;
01134   else
01135     return -a;
01136 }
01137 
01138 template <typename T>
01139 Taylor<T>
01140 abs(const Taylor<T>& a)
01141 {
01142   if (a.coeff(0) >= 0)
01143     return a;
01144   else
01145     return -a;
01146 }
01147 
01148 template <typename T>
01149 Taylor<T>
01150 max(const Taylor<T>& a,
01151     const Taylor<T>& b)
01152 {
01153   if (a.coeff(0) >= b.coeff(0))
01154     return a;
01155   else
01156     return b;
01157 }
01158 
01159 template <typename T>
01160 Taylor<T>
01161 max(const T& a,
01162     const Taylor<T>& b)
01163 {
01164   if (a >= b.coeff(0))
01165     return Taylor<T>(b.degree(), a);
01166   else
01167     return b;
01168 }
01169 
01170 template <typename T>
01171 Taylor<T>
01172 max(const Taylor<T>& a,
01173     const T& b)
01174 {
01175   if (a.coeff(0) >= b)
01176     return a;
01177   else
01178     return Taylor<T>(a.degree(), b);
01179 }
01180 
01181 template <typename T>
01182 Taylor<T>
01183 min(const Taylor<T>& a,
01184     const Taylor<T>& b)
01185 {
01186   if (a.coeff(0) <= b.coeff(0))
01187     return a;
01188   else
01189     return b;
01190 }
01191 
01192 template <typename T>
01193 Taylor<T>
01194 min(const T& a,
01195     const Taylor<T>& b)
01196 {
01197   if (a <= b.coeff(0))
01198     return Taylor<T>(b.degree(), a);
01199   else
01200     return b;
01201 }
01202 
01203 template <typename T>
01204 Taylor<T>
01205 min(const Taylor<T>& a,
01206     const T& b)
01207 {
01208   if (a.coeff(0) <= b)
01209     return a;
01210   else
01211     return Taylor<T>(a.degree(), b);
01212 }
01213 
01214 template <typename T>
01215 bool
01216 operator==(const Taylor<T>& a, 
01217      const Taylor<T>& b)
01218 {
01219   return a.coeff(0) == b.coeff(0);
01220 }
01221 
01222 template <typename T>
01223 bool
01224 operator==(const T& a, 
01225      const Taylor<T>& b)
01226 {
01227   return a == b.coeff(0);
01228 }
01229 
01230 template <typename T>
01231 bool
01232 operator==(const Taylor<T>& a, 
01233      const T& b)
01234 {
01235   return a.coeff(0) == b;
01236 }
01237 
01238 template <typename T>
01239 bool
01240 operator!=(const Taylor<T>& a, 
01241      const Taylor<T>& b)
01242 {
01243   return a.coeff(0) != b.coeff(0);
01244 }
01245 
01246 template <typename T>
01247 bool
01248 operator!=(const T& a, 
01249      const Taylor<T>& b)
01250 {
01251   return a != b.coeff(0);
01252 }
01253 
01254 template <typename T>
01255 bool
01256 operator!=(const Taylor<T>& a, 
01257      const T& b)
01258 {
01259   return a.coeff(0) != b;
01260 }
01261 
01262 template <typename T>
01263 bool
01264 operator<=(const Taylor<T>& a, 
01265      const Taylor<T>& b)
01266 {
01267   return a.coeff(0) <= b.coeff(0);
01268 }
01269 
01270 template <typename T>
01271 bool
01272 operator<=(const T& a, 
01273      const Taylor<T>& b)
01274 {
01275   return a <= b.coeff(0);
01276 }
01277 
01278 template <typename T>
01279 bool
01280 operator<=(const Taylor<T>& a, 
01281      const T& b)
01282 {
01283   return a.coeff(0) <= b;
01284 }
01285 
01286 template <typename T>
01287 bool
01288 operator>=(const Taylor<T>& a, 
01289      const Taylor<T>& b)
01290 {
01291   return a.coeff(0) >= b.coeff(0);
01292 }
01293 
01294 template <typename T>
01295 bool
01296 operator>=(const T& a, 
01297      const Taylor<T>& b)
01298 {
01299   return a >= b.coeff(0);
01300 }
01301 
01302 template <typename T>
01303 bool
01304 operator>=(const Taylor<T>& a, 
01305      const T& b)
01306 {
01307   return a.coeff(0) >= b;
01308 }
01309 
01310 template <typename T>
01311 bool
01312 operator<(const Taylor<T>& a, 
01313     const Taylor<T>& b)
01314 {
01315   return a.coeff(0) < b.coeff(0);
01316 }
01317 
01318 template <typename T>
01319 bool
01320 operator<(const T& a, 
01321     const Taylor<T>& b)
01322 {
01323   return a < b.coeff(0);
01324 }
01325 
01326 template <typename T>
01327 bool
01328 operator<(const Taylor<T>& a, 
01329     const T& b)
01330 {
01331   return a.coeff(0) < b;
01332 }
01333 
01334 template <typename T>
01335 bool
01336 operator>(const Taylor<T>& a, 
01337     const Taylor<T>& b)
01338 {
01339   return a.coeff(0) > b.coeff(0);
01340 }
01341 
01342 template <typename T>
01343 bool
01344 operator>(const T& a, 
01345     const Taylor<T>& b)
01346 {
01347   return a > b.coeff(0);
01348 }
01349 
01350 template <typename T>
01351 bool
01352 operator>(const Taylor<T>& a, 
01353     const T& b)
01354 {
01355   return a.coeff(0) > b;
01356 }
01357 
01358 template <typename T>
01359 bool toBool(const Taylor<T>& x) {
01360   bool is_zero = true;
01361   for (unsigned int i=0; i<=x.degree(); i++)
01362     is_zero = is_zero && (x.coeff(i) == 0.0);
01363   return !is_zero;
01364 }
01365 
01366 template <typename T>
01367 inline bool
01368 operator && (const Taylor<T>& x1, const Taylor<T>& x2)
01369 {
01370   return toBool(x1) && toBool(x2);
01371 }
01372 
01373 template <typename T>
01374 inline bool
01375 operator && (const typename Taylor<T>::value_type& a, const Taylor<T>& x2)
01376 {
01377   return a && toBool(x2);
01378 }
01379 
01380 template <typename T>
01381 inline bool
01382 operator && (const Taylor<T>& x1, const typename Taylor<T>::value_type& b)
01383 {
01384   return toBool(x1) && b;
01385 }
01386 
01387 template <typename T>
01388 inline bool
01389 operator || (const Taylor<T>& x1, const Taylor<T>& x2)
01390 {
01391   return toBool(x1) || toBool(x2);
01392 }
01393 
01394 template <typename T>
01395 inline bool
01396 operator || (const typename Taylor<T>::value_type& a, const Taylor<T>& x2)
01397 {
01398   return a || toBool(x2);
01399 }
01400 
01401 template <typename T>
01402 inline bool
01403 operator || (const Taylor<T>& x1, const typename Taylor<T>::value_type& b)
01404 {
01405   return toBool(x1) || b;
01406 }
01407 
01408 template <typename T>
01409 std::ostream& 
01410 operator << (std::ostream& os, const Taylor<T>& a)
01411 {
01412   os << "[ ";
01413       
01414   for (unsigned int i=0; i<=a.degree(); i++) {
01415     os << a.coeff(i) << " ";
01416   }
01417 
01418   os << "]";
01419   return os;
01420 }
01421 
01422 } // namespace Tay
01423 } // namespace Sacado
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines