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