Sacado_Tay_CacheTaylorOps.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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
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 #ifndef SACADO_TAY_CACHETAYLOROPS_HPP
00033 #define SACADO_TAY_CACHETAYLOROPS_HPP
00034 
00035 #include "Sacado_Tay_CacheTaylorExpr.hpp"
00036 
00037 #include <cmath>
00038 #include <valarray>
00039 #include <algorithm>  // for std::min and std::max
00040 #include <ostream>  // for std::ostream
00041 
00042 namespace Sacado {
00043 
00044   namespace Tay {
00045 
00046     // ---------------------- Unary Addition operator ------------------------
00047 
00048     template <typename ExprT>
00049     class UnaryPlusOp {
00050     public:
00051 
00052       typedef typename ExprT::value_type value_type;
00053 
00054       UnaryPlusOp(const ExprT& expr) {}
00055 
00056       void allocateCache(unsigned int d) const {}
00057 
00058       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00059   return expr.coeff(i);
00060       }
00061 
00062       value_type computeFastAccessCoeff(unsigned int i, 
00063           const ExprT& expr) const {
00064   return expr.fastAccessCoeff(i);
00065       }
00066 
00067     }; // class UnaryPlusOp
00068 
00069     // ---------------------- Unary Subtraction operator ---------------------
00070 
00071     template <typename ExprT>
00072     class UnaryMinusOp {
00073     public:
00074 
00075       typedef typename ExprT::value_type value_type;
00076 
00077       UnaryMinusOp(const ExprT& expr) {}
00078 
00079       void allocateCache(unsigned int d) const {}
00080 
00081       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00082   return -expr.coeff(i);
00083       }
00084 
00085       value_type computeFastAccessCoeff(unsigned int i, 
00086           const ExprT& expr) const {
00087   return -expr.fastAccessCoeff(i);
00088       }
00089 
00090     }; // class UnaryPlusOp
00091 
00092     // -------------------------- exp() function -----------------------------
00093 
00094     template <typename ExprT>
00095     class ExpOp {
00096     public:
00097 
00098       typedef typename ExprT::value_type value_type;
00099 
00100       ExpOp(const ExprT& expr) :
00101   c(),
00102   dc(-1) {}
00103 
00104       void allocateCache(unsigned int d) const { 
00105   c.resize(d+1,value_type(0));
00106       }
00107 
00108       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00109   if (static_cast<int>(i) > dc) {
00110     if (dc < 0) {
00111       c[0] = std::exp(expr.fastAccessCoeff(0));
00112       dc = 0;
00113     }
00114     for (unsigned int k=dc+1; k<=i; k++) {
00115       for (unsigned int j=1; j<=k; j++)
00116         c[k] += value_type(j)*c[k-j]*expr.coeff(j);
00117       c[k] /= value_type(k);
00118     }
00119     dc = i;
00120   }
00121   return c[i];
00122       }
00123 
00124       value_type computeFastAccessCoeff(unsigned int i, 
00125           const ExprT& expr) const
00126       {
00127   if (static_cast<int>(i) > dc) {
00128     if (dc < 0) {
00129       c[0] = std::exp(expr.fastAccessCoeff(0));
00130       dc = 0;
00131     }
00132     for (unsigned int k=dc+1; k<=i; k++) {
00133       for (unsigned int j=1; j<=k; j++)
00134         c[k] += value_type(j)*c[k-j]*expr.fastAccessCoeff(j);
00135       c[k] /= value_type(k);
00136     }
00137     dc = i;
00138   }
00139   return c[i];
00140       }
00141 
00142     protected:
00143 
00144       mutable std::valarray<value_type> c;
00145       mutable int dc;
00146 
00147     }; // class ExpOp
00148 
00149     // -------------------------- log() function -----------------------------
00150 
00151     template <typename ExprT>
00152     class LogOp {
00153     public:
00154 
00155       typedef typename ExprT::value_type value_type;
00156 
00157       LogOp(const ExprT& expr) :
00158   c(),
00159   dc(-1) 
00160       {}
00161 
00162       void allocateCache(unsigned int d) const { 
00163   c.resize(d+1,value_type(0));
00164       }
00165 
00166       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00167   if (static_cast<int>(i) > dc) {
00168     if (dc < 0) {
00169       c[0] = std::log(expr.fastAccessCoeff(0));
00170       dc = 0;
00171     }
00172     for (unsigned int k=dc+1; k<=i; k++) {
00173       c[k] = value_type(k)*expr.coeff(k);
00174       for (unsigned int j=1; j<=k-1; j++)
00175         c[k] -= value_type(j)*expr.coeff(k-j)*c[j];
00176       c[k] /= (value_type(k)*expr.fastAccessCoeff(0));
00177     }
00178     dc = i;
00179   }
00180   return c[i];
00181       }
00182 
00183       value_type computeFastAccessCoeff(unsigned int i, 
00184           const ExprT& expr) const
00185       {
00186   if (static_cast<int>(i) > dc) {
00187     if (dc < 0) {
00188       c[0] = std::log(expr.fastAccessCoeff(0));
00189       dc = 0;
00190     }
00191     for (unsigned int k=dc+1; k<=i; k++) {
00192       c[k] = value_type(k)*expr.fastAccessCoeff(k);
00193       for (unsigned int j=1; j<=k-1; j++)
00194         c[k] -= value_type(j)*expr.fastAccessCoeff(k-j)*c[j];
00195       c[k] /= (value_type(k)*expr.fastAccessCoeff(0));
00196     }
00197     dc = i;
00198   }
00199   return c[i];
00200       }
00201 
00202     protected:
00203 
00204       mutable std::valarray<value_type> c;
00205       mutable int dc;
00206 
00207     }; // class LogOp
00208 
00209     // -------------------------- sqrt() function -----------------------------
00210 
00211     template <typename ExprT>
00212     class SqrtOp {
00213     public:
00214 
00215       typedef typename ExprT::value_type value_type;
00216 
00217       SqrtOp(const ExprT& expr) :
00218   c(),
00219   dc(-1) {}
00220 
00221       void allocateCache(unsigned int d) const { 
00222   c.resize(d+1,value_type(0));
00223       }
00224 
00225       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00226   if (static_cast<int>(i) > dc) {
00227     if (dc < 0) {
00228       c[0] = std::sqrt(expr.fastAccessCoeff(0));
00229       dc = 0;
00230     }
00231     value_type tmp = value_type(2)*c[0];
00232     for (unsigned int k=dc+1; k<=i; k++) {
00233       c[k] = expr.coeff(k);
00234       for (unsigned int j=1; j<=k-1; j++)
00235         c[k] -= c[j]*c[k-j];
00236       c[k] /= tmp;
00237     }
00238     dc = i;
00239   }
00240   return c[i];
00241       }
00242 
00243       value_type computeFastAccessCoeff(unsigned int i, 
00244           const ExprT& expr) const
00245       {
00246   if (static_cast<int>(i) > dc) {
00247     if (dc < 0) {
00248       c[0] = std::sqrt(expr.fastAccessCoeff(0));
00249       dc = 0;
00250     }
00251     value_type tmp = value_type(2)*c[0];
00252     for (unsigned int k=dc+1; k<=i; k++) {
00253       c[k] = expr.fastAccessCoeff(k);
00254       for (unsigned int j=1; j<=k-1; j++)
00255         c[k] -= c[j]*c[k-j];
00256       c[k] /= tmp;
00257     }
00258     dc = i;
00259   }
00260   return c[i];
00261       }
00262 
00263     protected:
00264 
00265       mutable std::valarray<value_type> c;
00266       mutable int dc;
00267 
00268     }; // class SqrtOp
00269 
00270     // -------------------------- cos() function -----------------------------
00271 
00272     template <typename ExprT>
00273     class CosOp {
00274     public:
00275 
00276       typedef typename ExprT::value_type value_type;
00277 
00278       CosOp(const ExprT& expr) :
00279   c(),
00280   s(),
00281   dc(-1) {}
00282 
00283       void allocateCache(unsigned int d) const { 
00284   c.resize(d+1,value_type(0));
00285   s.resize(d+1,value_type(0));
00286       }
00287 
00288       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00289   if (static_cast<int>(i) > dc) {
00290     if (dc < 0) {
00291       c[0] = std::cos(expr.fastAccessCoeff(0));
00292       s[0] = std::sin(expr.fastAccessCoeff(0));
00293       dc = 0;
00294     }
00295     for (unsigned int k=dc+1; k<=i; k++) {
00296       for (unsigned int j=1; j<=k; j++) {
00297         c[k] -= value_type(j)*expr.coeff(j)*s[k-j];
00298         s[k] += value_type(j)*expr.coeff(j)*c[k-j];
00299       }
00300       c[k] /= value_type(k);
00301       s[k] /= value_type(k);
00302     }
00303     dc = i;
00304   }
00305   return c[i];
00306       }
00307 
00308       value_type computeFastAccessCoeff(unsigned int i, 
00309           const ExprT& expr) const
00310       {
00311   if (static_cast<int>(i) > dc) {
00312     if (dc < 0) {
00313       c[0] = std::cos(expr.fastAccessCoeff(0));
00314       s[0] = std::sin(expr.fastAccessCoeff(0));
00315       dc = 0;
00316     }
00317     for (unsigned int k=dc+1; k<=i; k++) {
00318       for (unsigned int j=1; j<=k; j++) {
00319         c[k] -= value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
00320         s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
00321       }
00322       c[k] /= value_type(k);
00323       s[k] /= value_type(k);
00324     }
00325     dc = i;
00326   }
00327   return c[i];
00328       }
00329 
00330     protected:
00331 
00332       mutable std::valarray<value_type> c;
00333       mutable std::valarray<value_type> s;
00334       mutable int dc;
00335 
00336     }; // class CosOp
00337 
00338     // -------------------------- sin() function -----------------------------
00339 
00340     template <typename ExprT>
00341     class SinOp {
00342     public:
00343 
00344       typedef typename ExprT::value_type value_type;
00345 
00346       SinOp(const ExprT& expr) :
00347   c(),
00348   s(),
00349   dc(-1) {}
00350 
00351       void allocateCache(unsigned int d) const { 
00352   c.resize(d+1,value_type(0));
00353   s.resize(d+1,value_type(0));
00354       }
00355 
00356       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00357   if (static_cast<int>(i) > dc) {
00358     if (dc < 0) {
00359       c[0] = std::cos(expr.fastAccessCoeff(0));
00360       s[0] = std::sin(expr.fastAccessCoeff(0));
00361       dc = 0;
00362     }
00363     for (unsigned int k=dc+1; k<=i; k++) {
00364       for (unsigned int j=1; j<=k; j++) {
00365         c[k] -= value_type(j)*expr.coeff(j)*s[k-j];
00366         s[k] += value_type(j)*expr.coeff(j)*c[k-j];
00367       }
00368       c[k] /= value_type(k);
00369       s[k] /= value_type(k);
00370     }
00371     dc = i;
00372   }
00373   return s[i];
00374       }
00375 
00376       value_type computeFastAccessCoeff(unsigned int i, 
00377           const ExprT& expr) const
00378       {
00379   if (static_cast<int>(i) > dc) {
00380     if (dc < 0) {
00381       c[0] = std::cos(expr.fastAccessCoeff(0));
00382       s[0] = std::sin(expr.fastAccessCoeff(0));
00383       dc = 0;
00384     }
00385     for (unsigned int k=dc+1; k<=i; k++) {
00386       for (unsigned int j=1; j<=k; j++) {
00387         c[k] -= value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
00388         s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
00389       }
00390       c[k] /= value_type(k);
00391       s[k] /= value_type(k);
00392     }
00393     dc = i;
00394   }
00395   return s[i];
00396       }
00397 
00398     protected:
00399 
00400       mutable std::valarray<value_type> c;
00401       mutable std::valarray<value_type> s;
00402       mutable int dc;
00403 
00404     }; // class SinOp
00405 
00406     // -------------------------- cosh() function -----------------------------
00407 
00408     template <typename ExprT>
00409     class CoshOp {
00410     public:
00411 
00412       typedef typename ExprT::value_type value_type;
00413 
00414       CoshOp(const ExprT& expr) :
00415   c(),
00416   s(),
00417   dc(-1) {}
00418 
00419       void allocateCache(unsigned int d) const { 
00420   c.resize(d+1,value_type(0));
00421   s.resize(d+1,value_type(0));
00422       }
00423 
00424       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00425   if (static_cast<int>(i) > dc) {
00426     if (dc < 0) {
00427       c[0] = std::cosh(expr.fastAccessCoeff(0));
00428       s[0] = std::sinh(expr.fastAccessCoeff(0));
00429       dc = 0;
00430     }
00431     for (unsigned int k=dc+1; k<=i; k++) {
00432       for (unsigned int j=1; j<=k; j++) {
00433         c[k] += value_type(j)*expr.coeff(j)*s[k-j];
00434         s[k] += value_type(j)*expr.coeff(j)*c[k-j];
00435       }
00436       c[k] /= value_type(k);
00437       s[k] /= value_type(k);
00438     }
00439     dc = i;
00440   }
00441   return c[i];
00442       }
00443 
00444       value_type computeFastAccessCoeff(unsigned int i, 
00445           const ExprT& expr) const
00446       {
00447   if (static_cast<int>(i) > dc) {
00448     if (dc < 0) {
00449       c[0] = std::cosh(expr.fastAccessCoeff(0));
00450       s[0] = std::sinh(expr.fastAccessCoeff(0));
00451       dc = 0;
00452     }
00453     for (unsigned int k=dc+1; k<=i; k++) {
00454       for (unsigned int j=1; j<=k; j++) {
00455         c[k] += value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
00456         s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
00457       }
00458       c[k] /= value_type(k);
00459       s[k] /= value_type(k);
00460     }
00461     dc = i;
00462   }
00463   return c[i];
00464       }
00465 
00466     protected:
00467 
00468       mutable std::valarray<value_type> c;
00469       mutable std::valarray<value_type> s;
00470       mutable int dc;
00471 
00472     }; // class CoshOp
00473 
00474     // -------------------------- sinh() function -----------------------------
00475 
00476     template <typename ExprT>
00477     class SinhOp {
00478     public:
00479 
00480       typedef typename ExprT::value_type value_type;
00481 
00482       SinhOp(const ExprT& expr) :
00483   c(),
00484   s(),
00485   dc(-1) {}
00486 
00487       void allocateCache(unsigned int d) const { 
00488   c.resize(d+1,value_type(0));
00489   s.resize(d+1,value_type(0));
00490       }
00491 
00492       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00493   if (static_cast<int>(i) > dc) {
00494     if (dc < 0) {
00495       c[0] = std::cosh(expr.fastAccessCoeff(0));
00496       s[0] = std::sinh(expr.fastAccessCoeff(0));
00497       dc = 0;
00498     }
00499     for (unsigned int k=dc+1; k<=i; k++) {
00500       for (unsigned int j=1; j<=k; j++) {
00501         c[k] += value_type(j)*expr.coeff(j)*s[k-j];
00502         s[k] += value_type(j)*expr.coeff(j)*c[k-j];
00503       }
00504       c[k] /= value_type(k);
00505       s[k] /= value_type(k);
00506     }
00507     dc = i;
00508   }
00509   return s[i];
00510       }
00511 
00512       value_type computeFastAccessCoeff(unsigned int i, 
00513           const ExprT& expr) const
00514       {
00515   if (static_cast<int>(i) > dc) {
00516     if (dc < 0) {
00517       c[0] = std::cosh(expr.fastAccessCoeff(0));
00518       s[0] = std::sinh(expr.fastAccessCoeff(0));
00519       dc = 0;
00520     }
00521     for (unsigned int k=dc+1; k<=i; k++) {
00522       for (unsigned int j=1; j<=k; j++) {
00523         c[k] += value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
00524         s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
00525       }
00526       c[k] /= value_type(k);
00527       s[k] /= value_type(k);
00528     }
00529     dc = i;
00530   }
00531   return s[i];
00532       }
00533 
00534     protected:
00535 
00536       mutable std::valarray<value_type> c;
00537       mutable std::valarray<value_type> s;
00538       mutable int dc;
00539 
00540     }; // class SinhOp
00541 
00542     // -------------------------- fabs() function -----------------------------
00543 
00544     template <typename ExprT>
00545     class FAbsOp {
00546     public:
00547 
00548       typedef typename ExprT::value_type value_type;
00549 
00550       FAbsOp(const ExprT& expr) {}
00551 
00552       void allocateCache(unsigned int d) const {}
00553 
00554       value_type computeCoeff(unsigned int i, const ExprT& expr) const {
00555   if (expr.fastAccessCoeff(0) > 0)
00556     return expr.coeff(i);
00557   else
00558     return -expr.coeff(i);
00559       }
00560 
00561       value_type computeFastAccessCoeff(unsigned int i, 
00562           const ExprT& expr) const
00563       {
00564   if (expr.fastAccessCoeff(0) > 0)
00565     return expr.fastAccessCoeff(i);
00566   else
00567     return -expr.fastAccessCoeff(i);
00568       }
00569 
00570     }; // class FAbsOp
00571 
00572   } // namespace Tay
00573 
00574 } // namespace Sacado
00575 
00576 #define TAYLOR_UNARYOP_MACRO(OPNAME,OP)         \
00577 namespace Sacado {              \
00578   namespace Tay {             \
00579     template <typename T>           \
00580     inline Expr< UnaryExpr< Expr<T>, OP > >       \
00581     OPNAME (const Expr<T>& expr)          \
00582     {                 \
00583       typedef UnaryExpr< Expr<T>, OP > expr_t;        \
00584                         \
00585       return Expr<expr_t>(expr_t(expr));        \
00586     }                 \
00587   }                 \
00588 }                                                                       \
00589                                                                         \
00590 namespace std {                                                         \
00591   using Sacado::Tay::OPNAME;                                            \
00592 }
00593 
00594 TAYLOR_UNARYOP_MACRO(operator+, UnaryPlusOp)
00595 TAYLOR_UNARYOP_MACRO(operator-, UnaryMinusOp)
00596 TAYLOR_UNARYOP_MACRO(exp, ExpOp)
00597 TAYLOR_UNARYOP_MACRO(log, LogOp)
00598 TAYLOR_UNARYOP_MACRO(sqrt, SqrtOp)
00599 TAYLOR_UNARYOP_MACRO(cos, CosOp)
00600 TAYLOR_UNARYOP_MACRO(sin, SinOp)
00601 TAYLOR_UNARYOP_MACRO(cosh, CoshOp)
00602 TAYLOR_UNARYOP_MACRO(sinh, SinhOp)
00603 TAYLOR_UNARYOP_MACRO(abs, FAbsOp)
00604 TAYLOR_UNARYOP_MACRO(fabs, FAbsOp)
00605 
00606 #undef TAYLOR_UNARYOP_MACRO
00607 
00608 namespace Sacado {
00609 
00610   namespace Tay {
00611 
00612     // ---------------------- Addition operator -----------------------------
00613 
00614     template <typename ExprT1, typename ExprT2>
00615     class AdditionOp {
00616     public:
00617       
00618       typedef typename ExprT1::value_type value_type_1;
00619       typedef typename ExprT2::value_type value_type_2;
00620       typedef typename Sacado::Promote<value_type_1,
00621                value_type_2>::type value_type;
00622     
00623       AdditionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00624 
00625       void allocateCache(unsigned int d) const {}
00626     
00627       value_type
00628       computeCoeff(unsigned int i, const ExprT1& expr1, 
00629        const ExprT2& expr2) const {
00630   return expr1.coeff(i) + expr2.coeff(i);
00631       }
00632       
00633       value_type
00634       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1, 
00635            const ExprT2& expr2) const {
00636   return expr1.fastAccessCoeff(i) + expr2.fastAccessCoeff(i);
00637       }
00638       
00639     }; // class AdditionOp
00640 
00641     template <typename ExprT1>
00642     class AdditionOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
00643     public:
00644 
00645       typedef typename ExprT1::value_type value_type;
00646       typedef ConstExpr<typename ExprT1::value_type> ExprT2;
00647 
00648       AdditionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00649 
00650       void allocateCache(unsigned int d) const {}
00651 
00652       value_type
00653       computeCoeff(unsigned int i, const ExprT1& expr1, 
00654        const ExprT2& expr2) const {
00655   if (i == 0)
00656     return expr1.coeff(i) + expr2.coeff(i);
00657   else
00658     return expr1.coeff(i);
00659       }
00660 
00661       value_type
00662       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00663            const ExprT2& expr2) const {
00664   if (i == 0)
00665     return expr1.fastAccessCoeff(i) + expr2.fastAccessCoeff(i);
00666   else
00667     return expr1.fastAccessCoeff(i);
00668       }
00669 
00670     }; // class AdditionOp
00671 
00672     template <typename ExprT2>
00673     class AdditionOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
00674     public:
00675 
00676       typedef typename ExprT2::value_type value_type;
00677       typedef ConstExpr<typename ExprT2::value_type> ExprT1;
00678 
00679       AdditionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00680 
00681       void allocateCache(unsigned int d) const {}
00682 
00683       value_type
00684       computeCoeff(unsigned int i, const ExprT1& expr1, 
00685        const ExprT2& expr2) const {
00686   if (i == 0)
00687     return expr1.coeff(i) + expr2.coeff(i);
00688   else
00689     return expr2.coeff(i);
00690       }
00691 
00692       value_type
00693       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00694            const ExprT2& expr2) const {
00695   if (i == 0)
00696     return expr1.fastAccessCoeff(i) + expr2.fastAccessCoeff(i);
00697   else
00698     return expr2.fastAccessCoeff(i);
00699       }
00700 
00701     }; // class AdditionOp
00702 
00703     // ---------------------- Subtraction operator ---------------------------
00704 
00705     template <typename ExprT1, typename ExprT2>
00706     class SubtractionOp {
00707     public:
00708       
00709       typedef typename ExprT1::value_type value_type_1;
00710       typedef typename ExprT2::value_type value_type_2;
00711       typedef typename Sacado::Promote<value_type_1,
00712                value_type_2>::type value_type;
00713     
00714       SubtractionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00715 
00716       void allocateCache(unsigned int d) const {}
00717     
00718       value_type
00719       computeCoeff(unsigned int i, const ExprT1& expr1, 
00720        const ExprT2& expr2) const {
00721   return expr1.coeff(i) - expr2.coeff(i);
00722       }
00723       
00724       value_type
00725       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1, 
00726            const ExprT2& expr2) const {
00727   return expr1.fastAccessCoeff(i) - expr2.fastAccessCoeff(i);
00728       }
00729       
00730     }; // class SubtractionOp
00731 
00732     template <typename ExprT1>
00733     class SubtractionOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
00734     public:
00735 
00736       typedef typename ExprT1::value_type value_type;
00737       typedef ConstExpr<typename ExprT1::value_type> ExprT2;
00738 
00739       SubtractionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00740 
00741       void allocateCache(unsigned int d) const {}
00742 
00743       value_type
00744       computeCoeff(unsigned int i, const ExprT1& expr1, 
00745        const ExprT2& expr2) const {
00746   if (i == 0)
00747     return expr1.coeff(i) - expr2.coeff(i);
00748   else
00749     return expr1.coeff(i);
00750       }
00751 
00752       value_type
00753       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00754            const ExprT2& expr2) const {
00755   if (i == 0)
00756     return expr1.fastAccessCoeff(i) - expr2.fastAccessCoeff(i);
00757   else
00758     return expr1.fastAccessCoeff(i);
00759       }
00760 
00761     }; // class SubtractionOp
00762 
00763     template <typename ExprT2>
00764     class SubtractionOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
00765     public:
00766 
00767       typedef typename ExprT2::value_type value_type;
00768       typedef ConstExpr<typename ExprT2::value_type> ExprT1;
00769 
00770       SubtractionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00771 
00772       void allocateCache(unsigned int d) const {}
00773 
00774       value_type
00775       computeCoeff(unsigned int i, const ExprT1& expr1, 
00776        const ExprT2& expr2) const {
00777   if (i == 0)
00778     return expr1.coeff(i) - expr2.coeff(i);
00779   else
00780     return -expr2.coeff(i);
00781       }
00782 
00783       value_type
00784       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00785            const ExprT2& expr2) const {
00786   if (i == 0)
00787     return expr1.fastAccessCoeff(i) - expr2.fastAccessCoeff(i);
00788   else
00789     return -expr2.fastAccessCoeff(i);
00790       }
00791 
00792     }; // class SubtractionOp
00793 
00794     // ---------------------- Multiplication operator -------------------------
00795 
00796     template <typename ExprT1, typename ExprT2>
00797     class MultiplicationOp {
00798     public:
00799       
00800       typedef typename ExprT1::value_type value_type_1;
00801       typedef typename ExprT2::value_type value_type_2;
00802       typedef typename Sacado::Promote<value_type_1,
00803                value_type_2>::type value_type;
00804     
00805       MultiplicationOp(const ExprT1& expr1, const ExprT2 expr2) :
00806   c(),
00807   dc(-1) {}
00808 
00809       void allocateCache(unsigned int d) const { 
00810   c.resize(d+1,value_type(0));
00811       }
00812     
00813       value_type
00814       computeCoeff(unsigned int i, const ExprT1& expr1, 
00815        const ExprT2& expr2) const {
00816   if (static_cast<int>(i) > dc) {
00817     for (unsigned int k=dc+1; k<=i; k++) {
00818       for (unsigned int j=0; j<=k; j++)
00819         c[k] += expr1.coeff(j)*expr2.coeff(k-j);
00820     }
00821     dc = i;
00822   }
00823   return c[i];
00824       }
00825       
00826       value_type
00827       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1, 
00828            const ExprT2& expr2) const {
00829   if (static_cast<int>(i) > dc) {
00830     for (unsigned int k=dc+1; k<=i; k++) {
00831       for (unsigned int j=0; j<=k; j++)
00832         c[k] += expr1.fastAccessCoeff(j)*expr2.fastAccessCoeff(k-j);
00833     }
00834     dc = i;
00835   }
00836   return c[i];
00837       }
00838 
00839     protected:
00840 
00841       mutable std::valarray<value_type> c;
00842       mutable int dc;
00843       
00844     }; // class MultiplicationOp
00845 
00846     template <typename ExprT1>
00847     class MultiplicationOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
00848     public:
00849 
00850       typedef typename ExprT1::value_type value_type;
00851       typedef ConstExpr<typename ExprT1::value_type> ExprT2;
00852 
00853       MultiplicationOp(const ExprT1& expr1, const ExprT2 expr2) {}
00854 
00855       void allocateCache(unsigned int d) const {}
00856 
00857       value_type
00858       computeCoeff(unsigned int i, const ExprT1& expr1, 
00859        const ExprT2& expr2) const {
00860   return expr1.coeff(i)*expr2.value();
00861       }
00862 
00863       value_type
00864       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00865            const ExprT2& expr2) const {
00866   return expr1.fastAccessCoeff(i)*expr2.value();
00867       }
00868 
00869     }; // class MultiplicationOp
00870 
00871     template <typename ExprT2>
00872     class MultiplicationOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
00873     public:
00874 
00875       typedef typename ExprT2::value_type value_type;
00876       typedef ConstExpr<typename ExprT2::value_type> ExprT1;
00877 
00878       MultiplicationOp(const ExprT1& expr1, const ExprT2 expr2) {}
00879 
00880       void allocateCache(unsigned int d) const {}
00881 
00882       value_type
00883       computeCoeff(unsigned int i, const ExprT1& expr1, 
00884        const ExprT2& expr2) const {
00885   return expr1.value()*expr2.coeff(i);
00886       }
00887 
00888       value_type
00889       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00890            const ExprT2& expr2) const {
00891   return expr1.value()*expr2.fastAccessCoeff(i);
00892       }
00893 
00894     }; // class MultiplicationOp
00895 
00896     // ---------------------- Division operator -------------------------
00897 
00898     template <typename ExprT1, typename ExprT2>
00899     class DivisionOp {
00900     public:
00901       
00902       typedef typename ExprT1::value_type value_type_1;
00903       typedef typename ExprT2::value_type value_type_2;
00904       typedef typename Sacado::Promote<value_type_1,
00905                value_type_2>::type value_type;
00906     
00907       DivisionOp(const ExprT1& expr1, const ExprT2 expr2) :
00908   c(),
00909   dc(-1) {}
00910 
00911       void allocateCache(unsigned int d) const { 
00912   c.resize(d+1,value_type(0));
00913       }
00914     
00915       value_type
00916       computeCoeff(unsigned int i, const ExprT1& expr1, 
00917        const ExprT2& expr2) const {
00918   if (static_cast<int>(i) > dc) {
00919     for (unsigned int k=dc+1; k<=i; k++) {
00920       c[k] = expr1.coeff(k);
00921       for (unsigned int j=1; j<=k; j++)
00922         c[k] -= expr2.coeff(j)*c[k-j];
00923       c[k] /= expr2.fastAccessCoeff(0);
00924     }
00925     dc = i;
00926   }
00927   return c[i];
00928       }
00929       
00930       value_type
00931       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1, 
00932            const ExprT2& expr2) const {
00933   if (static_cast<int>(i) > dc) {
00934     for (unsigned int k=dc+1; k<=i; k++) {
00935       c[k] = expr1.coeff(k);
00936       for (unsigned int j=1; j<=k; j++)
00937         c[k] -= expr2.fastAccessCoeff(j)*c[k-j];
00938       c[k] /= expr2.fastAccessCoeff(0);
00939     }
00940     dc = i;
00941   }
00942   return c[i];
00943       }
00944 
00945     protected:
00946 
00947       mutable std::valarray<value_type> c;
00948       mutable int dc;
00949       
00950     }; // class DivisionOp
00951 
00952     template <typename ExprT1>
00953     class DivisionOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
00954     public:
00955 
00956       typedef typename ExprT1::value_type value_type;
00957       typedef ConstExpr<typename ExprT1::value_type> ExprT2;
00958 
00959       DivisionOp(const ExprT1& expr1, const ExprT2 expr2) {}
00960 
00961       void allocateCache(unsigned int d) const {}
00962 
00963       value_type
00964       computeCoeff(unsigned int i, const ExprT1& expr1, 
00965        const ExprT2& expr2) const {
00966   return expr1.coeff(i)/expr2.value();
00967       }
00968 
00969       value_type
00970       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
00971            const ExprT2& expr2) const {
00972   return expr1.fastAccessCoeff(i)/expr2.value();
00973       }
00974 
00975     }; // class DivisionOp
00976 
00977     template <typename ExprT2>
00978     class DivisionOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
00979     public:
00980 
00981       typedef typename ExprT2::value_type value_type;
00982       typedef ConstExpr<typename ExprT2::value_type> ExprT1;
00983 
00984       DivisionOp(const ExprT1& expr1, const ExprT2 expr2) :
00985   c(),
00986   dc(-1) {}
00987 
00988       void allocateCache(unsigned int d) const { 
00989   c.resize(d+1,value_type(0));
00990       }
00991 
00992       value_type
00993       computeCoeff(unsigned int i, const ExprT1& expr1, 
00994        const ExprT2& expr2) const {
00995   if (static_cast<int>(i) > dc) {
00996     if (dc < 0) {
00997       c[0] = expr1.fastAccessCoeff(0) / expr2.fastAccessCoeff(0);
00998       dc = 0;
00999     }
01000     for (unsigned int k=dc+1; k<=i; k++) {
01001       for (unsigned int j=1; j<=k; j++)
01002         c[k] -= expr2.coeff(j)*c[k-j];
01003       c[k] /= expr2.fastAccessCoeff(0);
01004     }
01005     dc = i;
01006   }
01007   return c[i];
01008       }
01009 
01010       value_type
01011       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
01012            const ExprT2& expr2) const {
01013   if (static_cast<int>(i) > dc) {
01014     if (dc < 0) {
01015       c[0] = expr1.fastAccessCoeff(0) / expr2.fastAccessCoeff(0);
01016       dc = 0;
01017     }
01018     for (unsigned int k=dc+1; k<=i; k++) {
01019       c[k] = expr1.coeff(k);
01020       for (unsigned int j=1; j<=k; j++)
01021         c[k] -= expr2.fastAccessCoeff(j)*c[k-j];
01022       c[k] /= expr2.fastAccessCoeff(0);
01023     }
01024     dc = i;
01025   }
01026   return c[i];
01027       }
01028 
01029     protected:
01030 
01031       mutable std::valarray<value_type> c;
01032       mutable int dc;
01033 
01034     }; // class DivisionOp
01035 
01036     // ---------------------- Max operator -----------------------------
01037 
01038     template <typename ExprT1, typename ExprT2>
01039     class MaxOp {
01040     public:
01041       
01042       typedef typename ExprT1::value_type value_type_1;
01043       typedef typename ExprT2::value_type value_type_2;
01044       typedef typename Sacado::Promote<value_type_1,
01045                value_type_2>::type value_type;
01046     
01047       MaxOp(const ExprT1& expr1, const ExprT2 expr2) {}
01048 
01049       void allocateCache(unsigned int d) const {}
01050     
01051       value_type
01052       computeCoeff(unsigned int i, const ExprT1& expr1, 
01053        const ExprT2& expr2) const {
01054   if (i == 0)
01055     return std::max(expr1.coeff(0), expr2.coeff(0));
01056   else
01057     return expr1.coeff(0) >= expr2.coeff(0) ? expr1.coeff(i) : 
01058       expr2.coeff(i);
01059       }
01060       
01061       value_type
01062       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1, 
01063            const ExprT2& expr2) const {
01064   if (i == 0)
01065     return std::max(expr1.fastAccessCoeff(0), expr2.fastAccessCoeff(0));
01066   else
01067     return expr1.fastAccessCoeff(0) >= expr2.fastAccessCoeff(0) ? 
01068       expr1.fastAccessoeff(i) : expr2.fastAccessCoeff(i);
01069       }
01070       
01071     }; // class MaxOp
01072 
01073     template <typename ExprT1>
01074     class MaxOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
01075     public:
01076 
01077       typedef typename ExprT1::value_type value_type;
01078       typedef ConstExpr<typename ExprT1::value_type> ExprT2;
01079 
01080       MaxOp(const ExprT1& expr1, const ExprT2 expr2) {}
01081 
01082       void allocateCache(unsigned int d) const {}
01083 
01084       value_type
01085       computeCoeff(unsigned int i, const ExprT1& expr1, 
01086        const ExprT2& expr2) const {
01087   if (i == 0)
01088     return std::max(expr1.coeff(0), expr2.value());
01089   else
01090     return expr1.coeff(0) >= expr2.value() ? expr1.coeff(i) : 
01091       value_type(0);
01092       }
01093 
01094       value_type
01095       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
01096            const ExprT2& expr2) const {
01097   if (i == 0)
01098     return std::max(expr1.fastAccessCoeff(0), expr2.value());
01099   else
01100     return expr1.fastAccessCoeff(0) >= expr2.value() ? 
01101       expr1.fastAccessCoeff(i) : value_type(0);
01102       }
01103 
01104     }; // class MaxOp
01105 
01106     template <typename ExprT2>
01107     class MaxOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
01108     public:
01109 
01110       typedef typename ExprT2::value_type value_type;
01111       typedef ConstExpr<typename ExprT2::value_type> ExprT1;
01112 
01113       MaxOp(const ExprT1& expr1, const ExprT2 expr2) {}
01114 
01115       void allocateCache(unsigned int d) const {}
01116 
01117       value_type
01118       computeCoeff(unsigned int i, const ExprT1& expr1, 
01119        const ExprT2& expr2) const {
01120   if (i == 0)
01121     return std::max(expr1.value(), expr2.coeff(0));
01122   else
01123     return expr1.value() >= expr2.coeff(0) ? value_type(0) : 
01124       expr2.coeff(i);
01125       }
01126 
01127       value_type
01128       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
01129            const ExprT2& expr2) const {
01130   if (i == 0)
01131     return std::max(expr1.value(), expr2.fastAccessCoeff(0));
01132   else
01133     return expr1.value() >= expr2.fastAccessCoeff(0) ? value_type(0) : 
01134       expr2.fastAccessCoeff(i);
01135       }
01136 
01137     }; // class MaxOp
01138 
01139     // ---------------------- Min operator -----------------------------
01140 
01141     template <typename ExprT1, typename ExprT2>
01142     class MinOp {
01143     public:
01144       
01145       typedef typename ExprT1::value_type value_type_1;
01146       typedef typename ExprT2::value_type value_type_2;
01147       typedef typename Sacado::Promote<value_type_1,
01148                value_type_2>::type value_type;
01149     
01150       MinOp(const ExprT1& expr1, const ExprT2 expr2) {}
01151 
01152       void allocateCache(unsigned int d) const {}
01153     
01154       value_type
01155       computeCoeff(unsigned int i, const ExprT1& expr1, 
01156        const ExprT2& expr2) const {
01157   if (i == 0)
01158     return min(expr1.coeff(0), expr2.coeff(0));
01159   else
01160     return expr1.coeff(0) <= expr2.coeff(0) ? expr1.coeff(i) : 
01161       expr2.coeff(i);
01162       }
01163       
01164       value_type
01165       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1, 
01166            const ExprT2& expr2) const {
01167   if (i == 0)
01168     return min(expr1.fastAccessCoeff(0), expr2.fastAccessCoeff(0));
01169   else
01170     return expr1.fastAccessCoeff(0) <= expr2.fastAccessCoeff(0) ? 
01171       expr1.fastAccessCoeff(i) : expr2.fastAccessCoeff(i);
01172       }
01173       
01174     }; // class MinOp
01175 
01176     template <typename ExprT1>
01177     class MinOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
01178     public:
01179 
01180       typedef typename ExprT1::value_type value_type;
01181       typedef ConstExpr<typename ExprT1::value_type> ExprT2;
01182 
01183       MinOp(const ExprT1& expr1, const ExprT2 expr2) {}
01184 
01185       void allocateCache(unsigned int d) const {}
01186 
01187       value_type
01188       computeCoeff(unsigned int i, const ExprT1& expr1, 
01189        const ExprT2& expr2) const {
01190   if (i == 0)
01191     return std::min(expr1.coeff(0), expr2.value());
01192   else
01193     return expr1.coeff(0) <= expr2.value() ? expr1.coeff(i) : 
01194       value_type(0);
01195       }
01196 
01197       value_type
01198       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
01199            const ExprT2& expr2) const {
01200   if (i == 0)
01201     return std::min(expr1.fastAccessCoeff(0), expr2.value());
01202   else
01203     return expr1.fastAccessCoeff(0) <= expr2.value() ? 
01204       expr1.fastAccessCoeff(i) : value_type(0);
01205       }
01206 
01207     }; // class MinOp
01208 
01209     template <typename ExprT2>
01210     class MinOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
01211     public:
01212 
01213       typedef typename ExprT2::value_type value_type;
01214       typedef ConstExpr<typename ExprT2::value_type> ExprT1;
01215 
01216       MinOp(const ExprT1& expr1, const ExprT2 expr2) {}
01217 
01218       void allocateCache(unsigned int d) const {}
01219 
01220       value_type
01221       computeCoeff(unsigned int i, const ExprT1& expr1, 
01222        const ExprT2& expr2) const {
01223   if (i == 0)
01224     return std::min(expr1.value(), expr2.coeff(0));
01225   else
01226     return expr1.value() <= expr2.coeff(0) ? value_type(0) : 
01227       expr2.coeff(i);
01228       }
01229 
01230       value_type
01231       computeFastAccessCoeff(unsigned int i, const ExprT1& expr1,
01232            const ExprT2& expr2) const {
01233   if (i == 0)
01234     return std::min(expr1.value(), expr2.fastAccessCoeff(0));
01235   else
01236     return expr1.value() <= expr2.fastAccessCoeff(0) ? value_type(0) : 
01237       expr2.fastAccessCoeff(i);
01238       }
01239 
01240     }; // class MinOp
01241 
01242     // ------------------------ quadrature function ---------------------------
01243 
01244     template <typename ExprT1, typename ExprT2>
01245     class ASinQuadOp {
01246     public:
01247 
01248       typedef typename ExprT1::value_type value_type_1;
01249       typedef typename ExprT2::value_type value_type_2;
01250       typedef typename Sacado::Promote<value_type_1,
01251                value_type_2>::type value_type;
01252 
01253       ASinQuadOp(const ExprT1& expr1, const ExprT2& expr2) :
01254   c(),
01255   dc(-1)
01256       {}
01257 
01258       void allocateCache(unsigned int d) const { 
01259   c.resize(d+1,value_type(0));
01260       }
01261 
01262       value_type computeCoeff(unsigned int i, const ExprT1& expr1,
01263             const ExprT2& expr2) const {
01264   if (static_cast<int>(i) > dc) {
01265     if (dc < 0) {
01266       c[0] = std::asin(expr1.fastAccessCoeff(0));
01267       dc = 0;
01268     }
01269     for (unsigned int k=dc+1; k<=i; k++) {
01270       for (unsigned int j=1; j<=k; j++)
01271         c[k] += value_type(j)*expr2.coeff(k-j)*expr1.coeff(j);
01272       c[k] /= value_type(k);
01273     }
01274     dc = i;
01275   }
01276   return c[i];
01277       }
01278 
01279       value_type computeFastAccessCoeff(unsigned int i, 
01280           const ExprT1& expr1,
01281           const ExprT2& expr2) const
01282       {
01283   if (static_cast<int>(i) > dc) {
01284     if (dc < 0) {
01285       c[0] = std::asin(expr1.fastAccessCoeff(0));
01286       dc = 0;
01287     }
01288     for (unsigned int k=dc+1; k<=i; k++) {
01289       for (unsigned int j=1; j<=k; j++)
01290         c[k] += value_type(j)*expr2.fastAccessCoeff(k-j)*
01291     expr1.fastAccessCoeff(j);
01292       c[k] /= value_type(k);
01293     }
01294     dc = i;
01295   }
01296   return c[i];
01297       }
01298 
01299     protected:
01300 
01301       mutable std::valarray<value_type> c;
01302       mutable int dc;
01303 
01304     }; // class ASinQuadOp
01305 
01306     template <typename ExprT1, typename ExprT2>
01307     class ACosQuadOp {
01308     public:
01309 
01310       typedef typename ExprT1::value_type value_type_1;
01311       typedef typename ExprT2::value_type value_type_2;
01312       typedef typename Sacado::Promote<value_type_1,
01313                value_type_2>::type value_type;
01314 
01315       ACosQuadOp(const ExprT1& expr1, const ExprT2& expr2) :
01316   c(),
01317   dc(-1)
01318       {}
01319 
01320       void allocateCache(unsigned int d) const { 
01321   c.resize(d+1,value_type(0));
01322       }
01323 
01324       value_type computeCoeff(unsigned int i, const ExprT1& expr1,
01325             const ExprT2& expr2) const {
01326   if (static_cast<int>(i) > dc) {
01327     if (dc < 0) {
01328       c[0] = std::acos(expr1.fastAccessCoeff(0));
01329       dc = 0;
01330     }
01331     for (unsigned int k=dc+1; k<=i; k++) {
01332       for (unsigned int j=1; j<=k; j++)
01333         c[k] += value_type(j)*expr2.coeff(k-j)*expr1.coeff(j);
01334       c[k] /= value_type(k);
01335     }
01336     dc = i;
01337   }
01338   return c[i];
01339       }
01340 
01341       value_type computeFastAccessCoeff(unsigned int i, 
01342           const ExprT1& expr1,
01343           const ExprT2& expr2) const
01344       {
01345   if (static_cast<int>(i) > dc) {
01346     if (dc < 0) {
01347       c[0] = std::acos(expr1.fastAccessCoeff(0));
01348       dc = 0;
01349     }
01350     for (unsigned int k=dc+1; k<=i; k++) {
01351       for (unsigned int j=1; j<=k; j++)
01352         c[k] += value_type(j)*expr2.fastAccessCoeff(k-j)*
01353     expr1.fastAccessCoeff(j);
01354       c[k] /= value_type(k);
01355     }
01356     dc = i;
01357   }
01358   return c[i];
01359       }
01360 
01361     protected:
01362 
01363       mutable std::valarray<value_type> c;
01364       mutable int dc;
01365 
01366     }; // class ACosQuadOp
01367 
01368     template <typename ExprT1, typename ExprT2>
01369     class ATanQuadOp {
01370     public:
01371 
01372       typedef typename ExprT1::value_type value_type_1;
01373       typedef typename ExprT2::value_type value_type_2;
01374       typedef typename Sacado::Promote<value_type_1,
01375                value_type_2>::type value_type;
01376 
01377       ATanQuadOp(const ExprT1& expr1, const ExprT2& expr2) :
01378   c(),
01379   dc(-1)
01380       {}
01381 
01382       void allocateCache(unsigned int d) const { 
01383   c.resize(d+1,value_type(0));
01384       }
01385 
01386       value_type computeCoeff(unsigned int i, const ExprT1& expr1,
01387             const ExprT2& expr2) const {
01388   if (static_cast<int>(i) > dc) {
01389     if (dc < 0) {
01390       c[0] = std::atan(expr1.fastAccessCoeff(0));
01391       dc = 0;
01392     }
01393     for (unsigned int k=dc+1; k<=i; k++) {
01394       for (unsigned int j=1; j<=k; j++)
01395         c[k] += value_type(j)*expr2.coeff(k-j)*expr1.coeff(j);
01396       c[k] /= value_type(k);
01397     }
01398     dc = i;
01399   }
01400   return c[i];
01401       }
01402 
01403       value_type computeFastAccessCoeff(unsigned int i, 
01404           const ExprT1& expr1,
01405           const ExprT2& expr2) const
01406       {
01407   if (static_cast<int>(i) > dc) {
01408     if (dc < 0) {
01409       c[0] = std::atan(expr1.fastAccessCoeff(0));
01410       dc = 0;
01411     }
01412     for (unsigned int k=dc+1; k<=i; k++) {
01413       for (unsigned int j=1; j<=k; j++)
01414         c[k] += value_type(j)*expr2.fastAccessCoeff(k-j)*
01415     expr1.fastAccessCoeff(j);
01416       c[k] /= value_type(k);
01417     }
01418     dc = i;
01419   }
01420   return c[i];
01421       }
01422 
01423     protected:
01424 
01425       mutable std::valarray<value_type> c;
01426       mutable int dc;
01427 
01428     }; // class ATanQuadOp
01429 
01430   } // namespace Tay
01431 
01432 } // namespace Sacado
01433 
01434 #define TAYLOR_BINARYOP_MACRO(OPNAME,OP)        \
01435 namespace Sacado {              \
01436   namespace Tay {             \
01437     template <typename T1, typename T2>         \
01438     inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, OP > >     \
01439     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
01440     {                 \
01441       typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;    \
01442                       \
01443       return Expr<expr_t>(expr_t(expr1, expr2));      \
01444     }                 \
01445                   \
01446     template <typename T>           \
01447     inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>, \
01448            Expr<T>, OP > >        \
01449     OPNAME (const typename Expr<T>::value_type& c,      \
01450       const Expr<T>& expr)          \
01451     {                 \
01452       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
01453       typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;     \
01454                   \
01455       return Expr<expr_t>(expr_t(ConstT(c), expr));     \
01456     }                 \
01457                   \
01458     template <typename T>           \
01459     inline Expr< BinaryExpr< Expr<T>,         \
01460            ConstExpr<typename Expr<T>::value_type>, \
01461            OP > >         \
01462     OPNAME (const Expr<T>& expr,          \
01463       const typename Expr<T>::value_type& c)      \
01464     {                 \
01465       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
01466       typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;     \
01467                   \
01468       return Expr<expr_t>(expr_t(expr, ConstT(c)));     \
01469     }                 \
01470   }                 \
01471 }
01472 
01473 TAYLOR_BINARYOP_MACRO(operator+, AdditionOp)
01474 TAYLOR_BINARYOP_MACRO(operator-, SubtractionOp)
01475 TAYLOR_BINARYOP_MACRO(operator*, MultiplicationOp)
01476 TAYLOR_BINARYOP_MACRO(operator/, DivisionOp)
01477 
01478 #undef TAYLOR_BINARYOP_MACRO
01479 
01480   // The general definition of max/min works for Taylor variables too, except
01481   // we need to add a case when the argument types are different.  This 
01482   // can't conflict with the general definition, so we need to use
01483   // Substitution Failure Is Not An Error
01484 #include "Sacado_mpl_disable_if.hpp"
01485 #include "Sacado_mpl_is_same.hpp"
01486 
01487 #define TAYLOR_SFINAE_BINARYOP_MACRO(OPNAME,OP)       \
01488 namespace Sacado {              \
01489   namespace Tay {             \
01490     template <typename T1, typename T2>         \
01491     inline                                                              \
01492     typename                                                            \
01493     mpl::disable_if< mpl::is_same<T1,T2>,                               \
01494                      Expr<BinaryExpr<Expr<T1>, Expr<T2>, OP> > >::type  \
01495     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
01496     {                 \
01497       typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;    \
01498                       \
01499       return Expr<expr_t>(expr_t(expr1, expr2));      \
01500     }                 \
01501                   \
01502     template <typename T>           \
01503     inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>, \
01504            Expr<T>, OP > >        \
01505     OPNAME (const typename Expr<T>::value_type& c,      \
01506       const Expr<T>& expr)          \
01507     {                 \
01508       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
01509       typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;     \
01510                   \
01511       return Expr<expr_t>(expr_t(ConstT(c), expr));     \
01512     }                 \
01513                   \
01514     template <typename T>           \
01515     inline Expr< BinaryExpr< Expr<T>,         \
01516            ConstExpr<typename Expr<T>::value_type>, \
01517            OP > >         \
01518     OPNAME (const Expr<T>& expr,          \
01519       const typename Expr<T>::value_type& c)      \
01520     {                 \
01521       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
01522       typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;     \
01523                   \
01524       return Expr<expr_t>(expr_t(expr, ConstT(c)));     \
01525     }                 \
01526   }                 \
01527 }
01528 
01529 TAYLOR_SFINAE_BINARYOP_MACRO(max, MaxOp)
01530 TAYLOR_SFINAE_BINARYOP_MACRO(min, MinOp)
01531 
01532 #undef TAYLOR_SFINAE_BINARYOP_MACRO
01533 
01534 namespace std {
01535   using Sacado::Tay::min;
01536   using Sacado::Tay::max;
01537 }
01538 
01539 namespace Sacado {
01540 
01541   namespace Tay {
01542 
01543     template <typename T1, typename T2>
01544     inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, ASinQuadOp > >
01545     asin_quad (const Expr<T1>& expr1, const Expr<T2>& expr2)
01546     {
01547       typedef BinaryExpr< Expr<T1>, Expr<T2>, ASinQuadOp > expr_t;
01548  
01549       return Expr<expr_t>(expr_t(expr1, expr2));
01550     }
01551 
01552     template <typename T1, typename T2>
01553     inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, ACosQuadOp > >
01554     acos_quad (const Expr<T1>& expr1, const Expr<T2>& expr2)
01555     {
01556       typedef BinaryExpr< Expr<T1>, Expr<T2>, ACosQuadOp > expr_t;
01557  
01558       return Expr<expr_t>(expr_t(expr1, expr2));
01559     }
01560 
01561     template <typename T1, typename T2>
01562     inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, ATanQuadOp > >
01563     atan_quad (const Expr<T1>& expr1, const Expr<T2>& expr2)
01564     {
01565       typedef BinaryExpr< Expr<T1>, Expr<T2>, ATanQuadOp > expr_t;
01566  
01567       return Expr<expr_t>(expr_t(expr1, expr2));
01568     }
01569 
01570     template <typename ExprT1, typename ExprT2>
01571     struct PowExprType {
01572       typedef UnaryExpr< ExprT1, LogOp > T3;
01573       typedef BinaryExpr< ExprT2, Expr<T3>, MultiplicationOp > T4;
01574       typedef UnaryExpr< Expr<T4>, ExpOp > T5;
01575       
01576       typedef Expr<T5> expr_type;
01577     };
01578 
01579      template <typename ExprT2>
01580      struct PowExprType< typename ExprT2::value_type, ExprT2 > {
01581        typedef typename ExprT2::value_type T1;
01582       typedef BinaryExpr< ExprT2, ConstExpr<T1>, MultiplicationOp > T4;
01583       typedef UnaryExpr< Expr<T4>, ExpOp > T5;
01584       
01585       typedef Expr<T5> expr_type;
01586     };
01587 
01588     template <typename ExprT1>
01589     struct PowExprType<ExprT1,typename ExprT1::value_type> {
01590       typedef typename ExprT1::value_type T2;
01591       typedef UnaryExpr< ExprT1, LogOp > T3;
01592       typedef BinaryExpr< ConstExpr<T2>, Expr<T3>, MultiplicationOp > T4;
01593       typedef UnaryExpr< Expr<T4>, ExpOp > T5;
01594       
01595       typedef Expr<T5> expr_type;
01596     };
01597 
01598     template <typename T1, typename T2>
01599     inline typename PowExprType< Expr<T1>, Expr<T2> >::expr_type
01600     pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
01601     {
01602       // pow(x,y) = exp(y*log(x))
01603       return exp(expr2*log(expr1));
01604     }
01605 
01606     template <typename T>
01607     inline typename PowExprType< typename Expr<T>::value_type, Expr<T> >::expr_type
01608     pow (const typename Expr<T>::value_type& c, const Expr<T>& expr)
01609     {
01610       typedef ConstExpr<typename Expr<T>::value_type> ConstT;
01611 
01612       // pow(x,y) = exp(y*log(x))
01613       return exp(expr*std::log(c));
01614     }
01615 
01616     template <typename T>
01617     inline typename PowExprType< Expr<T>, typename Expr<T>::value_type >::expr_type
01618     pow (const Expr<T>& expr, const typename Expr<T>::value_type& c)
01619     {
01620       typedef ConstExpr<typename Expr<T>::value_type> ConstT;
01621 
01622       // pow(x,y) = exp(y*log(x))
01623       return exp(c*log(expr));
01624     }
01625 
01626     template <typename T>
01627     struct Log10ExprType {
01628       typedef UnaryExpr< Expr<T>, LogOp > T1;
01629       typedef ConstExpr<typename Expr<T>::value_type> ConstT;
01630       typedef BinaryExpr< Expr<T1>, ConstT, DivisionOp > T2;
01631       typedef Expr<T2> expr_type;
01632     };
01633 
01634     template <typename T>
01635     inline typename Log10ExprType<T>::expr_type
01636     log10 (const Expr<T>& expr)
01637     {
01638       // log10(x) = log(x)/log(10)
01639       return log(expr)/std::log(typename Expr<T>::value_type(10));
01640     }
01641 
01642     template <typename T>
01643     struct TanExprType {
01644       typedef UnaryExpr< Expr<T>, SinOp > T1;
01645       typedef UnaryExpr< Expr<T>, CosOp > T2;
01646       typedef BinaryExpr< Expr<T1>, Expr<T2>, DivisionOp > T3;
01647       typedef Expr<T3> expr_type;
01648     };
01649 
01650     template <typename T>
01651     inline typename TanExprType<T>::expr_type
01652     tan (const Expr<T>& expr)
01653     {
01654       // tan(x) = sin(x)/cos(x)
01655       return sin(expr)/cos(expr);
01656     }
01657 
01658     template <typename T>
01659     struct ASinExprType {
01660       typedef BinaryExpr< Expr<T>, Expr<T>, MultiplicationOp > T1;
01661       typedef ConstExpr<typename Expr<T>::value_type> ConstT;
01662       typedef BinaryExpr< ConstT, Expr<T1>, SubtractionOp > T2;
01663       typedef UnaryExpr< Expr<T2>, SqrtOp > T3;
01664       typedef BinaryExpr< ConstT, Expr<T3>, DivisionOp > T4;
01665       typedef BinaryExpr< Expr<T>, Expr<T4>, ASinQuadOp > T5;
01666       typedef Expr<T5> expr_type;
01667     };
01668 
01669     template <typename T>
01670     inline typename ASinExprType<T>::expr_type
01671     asin (const Expr<T>& expr)
01672     {
01673       typedef typename Expr<T>::value_type value_type;
01674 
01675       // asin(x) = integral of 1/sqrt(1-x*x)
01676       return asin_quad(expr, value_type(1)/sqrt(value_type(1)-expr*expr));
01677     }
01678 
01679     template <typename T>
01680     struct ACosExprType {
01681       typedef BinaryExpr< Expr<T>, Expr<T>, MultiplicationOp > T1;
01682       typedef ConstExpr<typename Expr<T>::value_type> ConstT;
01683       typedef BinaryExpr< ConstT, Expr<T1>, SubtractionOp > T2;
01684       typedef UnaryExpr< Expr<T2>, SqrtOp > T3;
01685       typedef BinaryExpr< ConstT, Expr<T3>, DivisionOp > T4;
01686       typedef BinaryExpr< Expr<T>, Expr<T4>, ACosQuadOp > T5;
01687       typedef Expr<T5> expr_type;
01688     };
01689 
01690     template <typename T>
01691     inline typename ACosExprType<T>::expr_type
01692     acos (const Expr<T>& expr)
01693     {
01694       typedef typename Expr<T>::value_type value_type;
01695 
01696       // acos(x) = integral of -1/sqrt(1-x*x)
01697       return acos_quad(expr, value_type(-1)/sqrt(value_type(1)-expr*expr));
01698     }
01699 
01700     template <typename T>
01701     struct ATanExprType {
01702       typedef BinaryExpr< Expr<T>, Expr<T>, MultiplicationOp > T1;
01703       typedef ConstExpr<typename Expr<T>::value_type> ConstT;
01704       typedef BinaryExpr< ConstT, Expr<T1>, AdditionOp > T2;
01705       typedef BinaryExpr< ConstT, Expr<T2>, DivisionOp > T3;
01706       typedef BinaryExpr< Expr<T>, Expr<T3>, ATanQuadOp > T4;
01707       typedef Expr<T4> expr_type;
01708     };
01709 
01710     template <typename T>
01711     inline typename ATanExprType<T>::expr_type
01712     atan (const Expr<T>& expr)
01713     {
01714       typedef typename Expr<T>::value_type value_type;
01715 
01716       // atan(x) = integral of 1/(1+x*x)
01717       return atan_quad(expr, value_type(1)/(value_type(1)+expr*expr));
01718     }
01719 
01720     template <typename T>
01721     struct TanhExprType {
01722       typedef UnaryExpr< Expr<T>, SinhOp > T1;
01723       typedef UnaryExpr< Expr<T>, CoshOp > T2;
01724       typedef BinaryExpr< Expr<T1>, Expr<T2>, DivisionOp > T3;
01725       typedef Expr<T3> expr_type;
01726     };
01727 
01728     template <typename T>
01729     inline typename TanhExprType<T>::expr_type
01730     tanh (const Expr<T>& expr)
01731     {
01732       // tanh(x) = sinh(x)/cosh(x)
01733       return sinh(expr)/cosh(expr);
01734     }
01735 
01736   } // namespace Tay
01737 
01738 } // namespace Sacado
01739 
01740 namespace std {
01741   using Sacado::Tay::pow;
01742   using Sacado::Tay::log10;
01743   using Sacado::Tay::tan;
01744   using Sacado::Tay::asin;
01745   using Sacado::Tay::acos;
01746   using Sacado::Tay::atan;
01747   using Sacado::Tay::tanh;
01748 }
01749 
01750 //-------------------------- Relational Operators -----------------------
01751 
01752 #define TAYLOR_RELOP_MACRO(OP)            \
01753 namespace Sacado {              \
01754   namespace Tay {             \
01755     template <typename ExprT1, typename ExprT2>       \
01756     inline bool               \
01757     operator OP (const Expr<ExprT1>& expr1,       \
01758      const Expr<ExprT2>& expr2)       \
01759     {                 \
01760       return expr1.fastAccessCoeff(0) OP expr2.fastAccessCoeff(0);  \
01761     }                 \
01762                   \
01763     template <typename ExprT2>            \
01764     inline bool               \
01765     operator OP (const typename Expr<ExprT2>::value_type& a,    \
01766      const Expr<ExprT2>& expr2)       \
01767     {                 \
01768       return a OP expr2.fastAccessCoeff(0);       \
01769     }                 \
01770                   \
01771     template <typename ExprT1>            \
01772     inline bool               \
01773     operator OP (const Expr<ExprT1>& expr1,       \
01774      const typename Expr<ExprT1>::value_type& b)    \
01775     {                 \
01776       return expr1.fastAccessCoeff(0) OP b;       \
01777     }                 \
01778   }                 \
01779 }
01780 
01781 TAYLOR_RELOP_MACRO(==)
01782 TAYLOR_RELOP_MACRO(!=)
01783 TAYLOR_RELOP_MACRO(<)
01784 TAYLOR_RELOP_MACRO(>)
01785 TAYLOR_RELOP_MACRO(<=)
01786 TAYLOR_RELOP_MACRO(>=)
01787 TAYLOR_RELOP_MACRO(<<=)
01788 TAYLOR_RELOP_MACRO(>>=)
01789 TAYLOR_RELOP_MACRO(&)
01790 TAYLOR_RELOP_MACRO(|)
01791 
01792 #undef TAYLOR_RELOP_MACRO
01793 
01794 namespace Sacado {
01795 
01796   namespace Tay {
01797 
01798     template <typename ExprT>
01799     inline bool operator ! (const Expr<ExprT>& expr) 
01800     {
01801       return ! expr.fastAccessCoeff(0);
01802     }
01803 
01804   } // namespace Tay
01805 
01806 } // namespace Sacado
01807 
01808 //-------------------------- I/O Operators -----------------------
01809 
01810 namespace Sacado {
01811 
01812   namespace Tay {
01813 
01814     template <typename ExprT>
01815     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
01816       os.setf(std::ios::fixed, std::ios::floatfield);
01817       os.width(12);
01818       os << "[";
01819       
01820       for (unsigned int i=0; i<=x.degree(); i++) {
01821         os.width(12);
01822         os << x.coeff(i);
01823       }
01824 
01825       os << "]";
01826       return os;
01827     }
01828 
01829   } // namespace Tay
01830 
01831 } // namespace Sacado
01832 
01833 
01834 #endif // SACADO_TAY_CACHETAYLOROPS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:19:36 2011 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.6.3