Sacado_Tay_CacheTaylorOps.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                           Sacado Package
//                 Copyright (2006) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
// (etphipp@sandia.gov).
//
// ***********************************************************************
// @HEADER

#ifndef SACADO_TAY_CACHETAYLOROPS_HPP
#define SACADO_TAY_CACHETAYLOROPS_HPP

#include "Sacado_Tay_CacheTaylorExpr.hpp"

#include <cmath>
#include <valarray>
#include <algorithm>    // for std::min and std::max
#include <ostream>      // for std::ostream

namespace Sacado {

  namespace Tay {

    // ---------------------- Unary Addition operator ------------------------

    template <typename ExprT>
    class UnaryPlusOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      UnaryPlusOp(const ExprT& expr) {}

      void allocateCache(int d) const {}

      value_type computeCoeff(int i, const ExprT& expr) const {
        return expr.coeff(i);
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const {
        return expr.fastAccessCoeff(i);
      }

    }; // class UnaryPlusOp

    // ---------------------- Unary Subtraction operator ---------------------

    template <typename ExprT>
    class UnaryMinusOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      UnaryMinusOp(const ExprT& expr) {}

      void allocateCache(int d) const {}

      value_type computeCoeff(int i, const ExprT& expr) const {
        return -expr.coeff(i);
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const {
        return -expr.fastAccessCoeff(i);
      }

    }; // class UnaryPlusOp

    // -------------------------- exp() function -----------------------------

    template <typename ExprT>
    class ExpOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      ExpOp(const ExprT& expr) :
        c(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::exp(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*c[k-j]*expr.coeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::exp(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*c[k-j]*expr.fastAccessCoeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class ExpOp

    // -------------------------- log() function -----------------------------

    template <typename ExprT>
    class LogOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      LogOp(const ExprT& expr) :
        c(),
        dc(-1)
      {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::log(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            c[k] = value_type(k)*expr.coeff(k);
            for (int j=1; j<=k-1; j++)
              c[k] -= value_type(j)*expr.coeff(k-j)*c[j];
            c[k] /= (value_type(k)*expr.fastAccessCoeff(0));
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::log(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            c[k] = value_type(k)*expr.fastAccessCoeff(k);
            for (int j=1; j<=k-1; j++)
              c[k] -= value_type(j)*expr.fastAccessCoeff(k-j)*c[j];
            c[k] /= (value_type(k)*expr.fastAccessCoeff(0));
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class LogOp

    // -------------------------- sqrt() function -----------------------------

    template <typename ExprT>
    class SqrtOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      SqrtOp(const ExprT& expr) :
        c(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::sqrt(expr.fastAccessCoeff(0));
            dc = 0;
          }
          value_type tmp = value_type(2)*c[0];
          for (int k=dc+1; k<=i; k++) {
            c[k] = expr.coeff(k);
            for (int j=1; j<=k-1; j++)
              c[k] -= c[j]*c[k-j];
            c[k] /= tmp;
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::sqrt(expr.fastAccessCoeff(0));
            dc = 0;
          }
          value_type tmp = value_type(2)*c[0];
          for (int k=dc+1; k<=i; k++) {
            c[k] = expr.fastAccessCoeff(k);
            for (int j=1; j<=k-1; j++)
              c[k] -= c[j]*c[k-j];
            c[k] /= tmp;
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class SqrtOp

    // -------------------------- cos() function -----------------------------

    template <typename ExprT>
    class CosOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      CosOp(const ExprT& expr) :
        c(),
        s(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
        s.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cos(expr.fastAccessCoeff(0));
            s[0] = std::sin(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] -= value_type(j)*expr.coeff(j)*s[k-j];
              s[k] += value_type(j)*expr.coeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cos(expr.fastAccessCoeff(0));
            s[0] = std::sin(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] -= value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
              s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable std::valarray<value_type> s;
      mutable int dc;

    }; // class CosOp

    // -------------------------- sin() function -----------------------------

    template <typename ExprT>
    class SinOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      SinOp(const ExprT& expr) :
        c(),
        s(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
        s.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cos(expr.fastAccessCoeff(0));
            s[0] = std::sin(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] -= value_type(j)*expr.coeff(j)*s[k-j];
              s[k] += value_type(j)*expr.coeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return s[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cos(expr.fastAccessCoeff(0));
            s[0] = std::sin(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] -= value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
              s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return s[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable std::valarray<value_type> s;
      mutable int dc;

    }; // class SinOp

    // -------------------------- cosh() function -----------------------------

    template <typename ExprT>
    class CoshOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      CoshOp(const ExprT& expr) :
        c(),
        s(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
        s.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cosh(expr.fastAccessCoeff(0));
            s[0] = std::sinh(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] += value_type(j)*expr.coeff(j)*s[k-j];
              s[k] += value_type(j)*expr.coeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cosh(expr.fastAccessCoeff(0));
            s[0] = std::sinh(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] += value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
              s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable std::valarray<value_type> s;
      mutable int dc;

    }; // class CoshOp

    // -------------------------- sinh() function -----------------------------

    template <typename ExprT>
    class SinhOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      SinhOp(const ExprT& expr) :
        c(),
        s(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
        s.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cosh(expr.fastAccessCoeff(0));
            s[0] = std::sinh(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] += value_type(j)*expr.coeff(j)*s[k-j];
              s[k] += value_type(j)*expr.coeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return s[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::cosh(expr.fastAccessCoeff(0));
            s[0] = std::sinh(expr.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++) {
              c[k] += value_type(j)*expr.fastAccessCoeff(j)*s[k-j];
              s[k] += value_type(j)*expr.fastAccessCoeff(j)*c[k-j];
            }
            c[k] /= value_type(k);
            s[k] /= value_type(k);
          }
          dc = i;
        }
        return s[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable std::valarray<value_type> s;
      mutable int dc;

    }; // class SinhOp

    // -------------------------- fabs() function -----------------------------

    template <typename ExprT>
    class FAbsOp {
    public:

      typedef typename ExprT::value_type value_type;
      typedef typename ExprT::scalar_type scalar_type;
      typedef typename ExprT::base_expr_type base_expr_type;

      FAbsOp(const ExprT& expr) {}

      void allocateCache(int d) const {}

      value_type computeCoeff(int i, const ExprT& expr) const {
        if (expr.fastAccessCoeff(0) > 0)
          return expr.coeff(i);
        else
          return -expr.coeff(i);
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT& expr) const
      {
        if (expr.fastAccessCoeff(0) > 0)
          return expr.fastAccessCoeff(i);
        else
          return -expr.fastAccessCoeff(i);
      }

    }; // class FAbsOp

  } // namespace Tay

} // namespace Sacado

#define TAYLOR_UNARYOP_MACRO(OPNAME,OP)                                 \
namespace Sacado {                                                      \
  namespace Tay {                                                       \
    template <typename T>                                               \
    inline Expr< UnaryExpr< Expr<T>, OP > >                             \
    OPNAME (const Expr<T>& expr)                                        \
    {                                                                   \
      typedef UnaryExpr< Expr<T>, OP > expr_t;                          \
                                                                        \
      return Expr<expr_t>(expr_t(expr));                                \
    }                                                                   \
  }                                                                     \
}                                                                       \
                                                                        \
namespace std {                                                         \
  using Sacado::Tay::OPNAME;                                            \
}

TAYLOR_UNARYOP_MACRO(operator+, UnaryPlusOp)
TAYLOR_UNARYOP_MACRO(operator-, UnaryMinusOp)
TAYLOR_UNARYOP_MACRO(exp, ExpOp)
TAYLOR_UNARYOP_MACRO(log, LogOp)
TAYLOR_UNARYOP_MACRO(sqrt, SqrtOp)
TAYLOR_UNARYOP_MACRO(cos, CosOp)
TAYLOR_UNARYOP_MACRO(sin, SinOp)
TAYLOR_UNARYOP_MACRO(cosh, CoshOp)
TAYLOR_UNARYOP_MACRO(sinh, SinhOp)
TAYLOR_UNARYOP_MACRO(abs, FAbsOp)
TAYLOR_UNARYOP_MACRO(fabs, FAbsOp)

#undef TAYLOR_UNARYOP_MACRO

namespace Sacado {

  namespace Tay {

    // ---------------------- Addition operator -----------------------------

    template <typename ExprT1, typename ExprT2>
    class AdditionOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      AdditionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        return expr1.coeff(i) + expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        return expr1.fastAccessCoeff(i) + expr2.fastAccessCoeff(i);
      }

    }; // class AdditionOp

    template <typename ExprT1>
    class AdditionOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
    public:

      typedef typename ExprT1::value_type value_type;
      typedef typename ExprT1::scalar_type scalar_type;
      typedef typename ExprT1::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT1::value_type> ExprT2;

      AdditionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return expr1.coeff(i) + expr2.coeff(i);
        else
          return expr1.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return expr1.fastAccessCoeff(i) + expr2.fastAccessCoeff(i);
        else
          return expr1.fastAccessCoeff(i);
      }

    }; // class AdditionOp

    template <typename ExprT2>
    class AdditionOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
    public:

      typedef typename ExprT2::value_type value_type;
      typedef typename ExprT2::scalar_type scalar_type;
      typedef typename ExprT2::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT2::value_type> ExprT1;

      AdditionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return expr1.coeff(i) + expr2.coeff(i);
        else
          return expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return expr1.fastAccessCoeff(i) + expr2.fastAccessCoeff(i);
        else
          return expr2.fastAccessCoeff(i);
      }

    }; // class AdditionOp

    // ---------------------- Subtraction operator ---------------------------

    template <typename ExprT1, typename ExprT2>
    class SubtractionOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      SubtractionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        return expr1.coeff(i) - expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        return expr1.fastAccessCoeff(i) - expr2.fastAccessCoeff(i);
      }

    }; // class SubtractionOp

    template <typename ExprT1>
    class SubtractionOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
    public:

      typedef typename ExprT1::value_type value_type;
      typedef typename ExprT1::scalar_type scalar_type;
      typedef typename ExprT1::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT1::value_type> ExprT2;

      SubtractionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return expr1.coeff(i) - expr2.coeff(i);
        else
          return expr1.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return expr1.fastAccessCoeff(i) - expr2.fastAccessCoeff(i);
        else
          return expr1.fastAccessCoeff(i);
      }

    }; // class SubtractionOp

    template <typename ExprT2>
    class SubtractionOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
    public:

      typedef typename ExprT2::value_type value_type;
      typedef typename ExprT2::scalar_type scalar_type;
      typedef typename ExprT2::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT2::value_type> ExprT1;

      SubtractionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return expr1.coeff(i) - expr2.coeff(i);
        else
          return -expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return expr1.fastAccessCoeff(i) - expr2.fastAccessCoeff(i);
        else
          return -expr2.fastAccessCoeff(i);
      }

    }; // class SubtractionOp

    // ---------------------- Multiplication operator -------------------------

    template <typename ExprT1, typename ExprT2>
    class MultiplicationOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      MultiplicationOp(const ExprT1& expr1, const ExprT2 expr2) :
        c(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          for (int k=dc+1; k<=i; k++) {
            for (int j=0; j<=k; j++)
              c[k] += expr1.coeff(j)*expr2.coeff(k-j);
          }
          dc = i;
        }
        return c[i];
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          for (int k=dc+1; k<=i; k++) {
            for (int j=0; j<=k; j++)
              c[k] += expr1.fastAccessCoeff(j)*expr2.fastAccessCoeff(k-j);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class MultiplicationOp

    template <typename ExprT1>
    class MultiplicationOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
    public:

      typedef typename ExprT1::value_type value_type;
      typedef typename ExprT1::scalar_type scalar_type;
      typedef typename ExprT1::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT1::value_type> ExprT2;

      MultiplicationOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        return expr1.coeff(i)*expr2.value();
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        return expr1.fastAccessCoeff(i)*expr2.value();
      }

    }; // class MultiplicationOp

    template <typename ExprT2>
    class MultiplicationOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
    public:

      typedef typename ExprT2::value_type value_type;
      typedef typename ExprT2::scalar_type scalar_type;
      typedef typename ExprT2::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT2::value_type> ExprT1;

      MultiplicationOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        return expr1.value()*expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        return expr1.value()*expr2.fastAccessCoeff(i);
      }

    }; // class MultiplicationOp

    // ---------------------- Division operator -------------------------

    template <typename ExprT1, typename ExprT2>
    class DivisionOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      DivisionOp(const ExprT1& expr1, const ExprT2 expr2) :
        c(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          for (int k=dc+1; k<=i; k++) {
            c[k] = expr1.coeff(k);
            for (int j=1; j<=k; j++)
              c[k] -= expr2.coeff(j)*c[k-j];
            c[k] /= expr2.fastAccessCoeff(0);
          }
          dc = i;
        }
        return c[i];
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          for (int k=dc+1; k<=i; k++) {
            c[k] = expr1.coeff(k);
            for (int j=1; j<=k; j++)
              c[k] -= expr2.fastAccessCoeff(j)*c[k-j];
            c[k] /= expr2.fastAccessCoeff(0);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class DivisionOp

    template <typename ExprT1>
    class DivisionOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
    public:

      typedef typename ExprT1::value_type value_type;
      typedef typename ExprT1::scalar_type scalar_type;
      typedef typename ExprT1::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT1::value_type> ExprT2;

      DivisionOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        return expr1.coeff(i)/expr2.value();
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        return expr1.fastAccessCoeff(i)/expr2.value();
      }

    }; // class DivisionOp

    template <typename ExprT2>
    class DivisionOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
    public:

      typedef typename ExprT2::value_type value_type;
      typedef typename ExprT2::scalar_type scalar_type;
      typedef typename ExprT2::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT2::value_type> ExprT1;

      DivisionOp(const ExprT1& expr1, const ExprT2 expr2) :
        c(),
        dc(-1) {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = expr1.fastAccessCoeff(0) / expr2.fastAccessCoeff(0);
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] -= expr2.coeff(j)*c[k-j];
            c[k] /= expr2.fastAccessCoeff(0);
          }
          dc = i;
        }
        return c[i];
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = expr1.fastAccessCoeff(0) / expr2.fastAccessCoeff(0);
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            c[k] = expr1.coeff(k);
            for (int j=1; j<=k; j++)
              c[k] -= expr2.fastAccessCoeff(j)*c[k-j];
            c[k] /= expr2.fastAccessCoeff(0);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class DivisionOp

    // ---------------------- Max operator -----------------------------

    template <typename ExprT1, typename ExprT2>
    class MaxOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      MaxOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return std::max(expr1.coeff(0), expr2.coeff(0));
        else
          return expr1.coeff(0) >= expr2.coeff(0) ? expr1.coeff(i) :
            expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return std::max(expr1.fastAccessCoeff(0), expr2.fastAccessCoeff(0));
        else
          return expr1.fastAccessCoeff(0) >= expr2.fastAccessCoeff(0) ?
            expr1.fastAccessoeff(i) : expr2.fastAccessCoeff(i);
      }

    }; // class MaxOp

    template <typename ExprT1>
    class MaxOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
    public:

      typedef typename ExprT1::value_type value_type;
      typedef typename ExprT1::scalar_type scalar_type;
      typedef typename ExprT1::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT1::value_type> ExprT2;

      MaxOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return std::max(expr1.coeff(0), expr2.value());
        else
          return expr1.coeff(0) >= expr2.value() ? expr1.coeff(i) :
            value_type(0);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return std::max(expr1.fastAccessCoeff(0), expr2.value());
        else
          return expr1.fastAccessCoeff(0) >= expr2.value() ?
            expr1.fastAccessCoeff(i) : value_type(0);
      }

    }; // class MaxOp

    template <typename ExprT2>
    class MaxOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
    public:

      typedef typename ExprT2::value_type value_type;
      typedef typename ExprT2::scalar_type scalar_type;
      typedef typename ExprT2::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT2::value_type> ExprT1;

      MaxOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return std::max(expr1.value(), expr2.coeff(0));
        else
          return expr1.value() >= expr2.coeff(0) ? value_type(0) :
            expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return std::max(expr1.value(), expr2.fastAccessCoeff(0));
        else
          return expr1.value() >= expr2.fastAccessCoeff(0) ? value_type(0) :
            expr2.fastAccessCoeff(i);
      }

    }; // class MaxOp

    // ---------------------- Min operator -----------------------------

    template <typename ExprT1, typename ExprT2>
    class MinOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;
      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      MinOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return min(expr1.coeff(0), expr2.coeff(0));
        else
          return expr1.coeff(0) <= expr2.coeff(0) ? expr1.coeff(i) :
            expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return min(expr1.fastAccessCoeff(0), expr2.fastAccessCoeff(0));
        else
          return expr1.fastAccessCoeff(0) <= expr2.fastAccessCoeff(0) ?
            expr1.fastAccessCoeff(i) : expr2.fastAccessCoeff(i);
      }

    }; // class MinOp

    template <typename ExprT1>
    class MinOp<ExprT1, ConstExpr<typename ExprT1::value_type> > {
    public:

      typedef typename ExprT1::value_type value_type;
      typedef typename ExprT1::scalar_type scalar_type;
      typedef typename ExprT1::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT1::value_type> ExprT2;

      MinOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return std::min(expr1.coeff(0), expr2.value());
        else
          return expr1.coeff(0) <= expr2.value() ? expr1.coeff(i) :
            value_type(0);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return std::min(expr1.fastAccessCoeff(0), expr2.value());
        else
          return expr1.fastAccessCoeff(0) <= expr2.value() ?
            expr1.fastAccessCoeff(i) : value_type(0);
      }

    }; // class MinOp

    template <typename ExprT2>
    class MinOp<ConstExpr<typename ExprT2::value_type>, ExprT2 > {
    public:

      typedef typename ExprT2::value_type value_type;
      typedef typename ExprT2::scalar_type scalar_type;
      typedef typename ExprT2::base_expr_type base_expr_type;
      typedef ConstExpr<typename ExprT2::value_type> ExprT1;

      MinOp(const ExprT1& expr1, const ExprT2 expr2) {}

      void allocateCache(int d) const {}

      value_type
      computeCoeff(int i, const ExprT1& expr1,
                   const ExprT2& expr2) const {
        if (i == 0)
          return std::min(expr1.value(), expr2.coeff(0));
        else
          return expr1.value() <= expr2.coeff(0) ? value_type(0) :
            expr2.coeff(i);
      }

      value_type
      computeFastAccessCoeff(int i, const ExprT1& expr1,
                             const ExprT2& expr2) const {
        if (i == 0)
          return std::min(expr1.value(), expr2.fastAccessCoeff(0));
        else
          return expr1.value() <= expr2.fastAccessCoeff(0) ? value_type(0) :
            expr2.fastAccessCoeff(i);
      }

    }; // class MinOp

    // ------------------------ quadrature function ---------------------------

    template <typename ExprT1, typename ExprT2>
    class ASinQuadOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      ASinQuadOp(const ExprT1& expr1, const ExprT2& expr2) :
        c(),
        dc(-1)
      {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT1& expr1,
                              const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::asin(expr1.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*expr2.coeff(k-j)*expr1.coeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT1& expr1,
                                        const ExprT2& expr2) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::asin(expr1.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*expr2.fastAccessCoeff(k-j)*
                expr1.fastAccessCoeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class ASinQuadOp

    template <typename ExprT1, typename ExprT2>
    class ACosQuadOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      ACosQuadOp(const ExprT1& expr1, const ExprT2& expr2) :
        c(),
        dc(-1)
      {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT1& expr1,
                              const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::acos(expr1.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*expr2.coeff(k-j)*expr1.coeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT1& expr1,
                                        const ExprT2& expr2) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::acos(expr1.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*expr2.fastAccessCoeff(k-j)*
                expr1.fastAccessCoeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class ACosQuadOp

    template <typename ExprT1, typename ExprT2>
    class ATanQuadOp {
    public:

      typedef typename ExprT1::value_type value_type_1;
      typedef typename ExprT2::value_type value_type_2;
      typedef typename Sacado::Promote<value_type_1,
                                       value_type_2>::type value_type;

      typedef typename ExprT1::scalar_type scalar_type_1;
      typedef typename ExprT2::scalar_type scalar_type_2;
      typedef typename Sacado::Promote<scalar_type_1,
                                       scalar_type_2>::type scalar_type;

      typedef typename ExprT1::base_expr_type base_expr_type_1;
      typedef typename ExprT2::base_expr_type base_expr_type_2;
      typedef typename Sacado::Promote<base_expr_type_1,
                                       base_expr_type_2>::type base_expr_type;

      ATanQuadOp(const ExprT1& expr1, const ExprT2& expr2) :
        c(),
        dc(-1)
      {}

      void allocateCache(int d) const {
        c.resize(d+1,value_type(0));
      }

      value_type computeCoeff(int i, const ExprT1& expr1,
                              const ExprT2& expr2) const {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::atan(expr1.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*expr2.coeff(k-j)*expr1.coeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

      value_type computeFastAccessCoeff(int i,
                                        const ExprT1& expr1,
                                        const ExprT2& expr2) const
      {
        if (static_cast<int>(i) > dc) {
          if (dc < 0) {
            c[0] = std::atan(expr1.fastAccessCoeff(0));
            dc = 0;
          }
          for (int k=dc+1; k<=i; k++) {
            for (int j=1; j<=k; j++)
              c[k] += value_type(j)*expr2.fastAccessCoeff(k-j)*
                expr1.fastAccessCoeff(j);
            c[k] /= value_type(k);
          }
          dc = i;
        }
        return c[i];
      }

    protected:

      mutable std::valarray<value_type> c;
      mutable int dc;

    }; // class ATanQuadOp

  } // namespace Tay

} // namespace Sacado

#define TAYLOR_BINARYOP_MACRO(OPNAME,OP)                                \
namespace Sacado {                                                      \
  namespace Tay {                                                       \
    template <typename T1, typename T2>                                 \
    inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, OP > >                 \
    OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)               \
    {                                                                   \
      typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;              \
                                                                        \
      return Expr<expr_t>(expr_t(expr1, expr2));                        \
    }                                                                   \
                                                                        \
    template <typename T>                                               \
    inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>,   \
                             Expr<T>, OP > >                            \
    OPNAME (const typename Expr<T>::value_type& c,                      \
            const Expr<T>& expr)                                        \
    {                                                                   \
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;           \
      typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;                 \
                                                                        \
      return Expr<expr_t>(expr_t(ConstT(c), expr));                     \
    }                                                                   \
                                                                        \
    template <typename T>                                               \
    inline Expr< BinaryExpr< Expr<T>,                                   \
                             ConstExpr<typename Expr<T>::value_type>,   \
                             OP > >                                     \
    OPNAME (const Expr<T>& expr,                                        \
            const typename Expr<T>::value_type& c)                      \
    {                                                                   \
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;           \
      typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;                 \
                                                                        \
      return Expr<expr_t>(expr_t(expr, ConstT(c)));                     \
    }                                                                   \
  }                                                                     \
}

TAYLOR_BINARYOP_MACRO(operator+, AdditionOp)
TAYLOR_BINARYOP_MACRO(operator-, SubtractionOp)
TAYLOR_BINARYOP_MACRO(operator*, MultiplicationOp)
TAYLOR_BINARYOP_MACRO(operator/, DivisionOp)

#undef TAYLOR_BINARYOP_MACRO

  // The general definition of max/min works for Taylor variables too, except
  // we need to add a case when the argument types are different.  This
  // can't conflict with the general definition, so we need to use
  // Substitution Failure Is Not An Error
#include "Sacado_mpl_disable_if.hpp"
#include "Sacado_mpl_is_same.hpp"

#define TAYLOR_SFINAE_BINARYOP_MACRO(OPNAME,OP)                         \
namespace Sacado {                                                      \
  namespace Tay {                                                       \
    template <typename T1, typename T2>                                 \
    inline                                                              \
    typename                                                            \
    mpl::disable_if< mpl::is_same<T1,T2>,                               \
                     Expr<BinaryExpr<Expr<T1>, Expr<T2>, OP> > >::type  \
    OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)               \
    {                                                                   \
      typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;              \
                                                                        \
      return Expr<expr_t>(expr_t(expr1, expr2));                        \
    }                                                                   \
                                                                        \
    template <typename T>                                               \
    inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>,   \
                             Expr<T>, OP > >                            \
    OPNAME (const typename Expr<T>::value_type& c,                      \
            const Expr<T>& expr)                                        \
    {                                                                   \
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;           \
      typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;                 \
                                                                        \
      return Expr<expr_t>(expr_t(ConstT(c), expr));                     \
    }                                                                   \
                                                                        \
    template <typename T>                                               \
    inline Expr< BinaryExpr< Expr<T>,                                   \
                             ConstExpr<typename Expr<T>::value_type>,   \
                             OP > >                                     \
    OPNAME (const Expr<T>& expr,                                        \
            const typename Expr<T>::value_type& c)                      \
    {                                                                   \
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;           \
      typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;                 \
                                                                        \
      return Expr<expr_t>(expr_t(expr, ConstT(c)));                     \
    }                                                                   \
  }                                                                     \
}

TAYLOR_SFINAE_BINARYOP_MACRO(max, MaxOp)
TAYLOR_SFINAE_BINARYOP_MACRO(min, MinOp)

#undef TAYLOR_SFINAE_BINARYOP_MACRO

namespace std {
  using Sacado::Tay::min;
  using Sacado::Tay::max;
}

namespace Sacado {

  namespace Tay {

    template <typename T1, typename T2>
    inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, ASinQuadOp > >
    asin_quad (const Expr<T1>& expr1, const Expr<T2>& expr2)
    {
      typedef BinaryExpr< Expr<T1>, Expr<T2>, ASinQuadOp > expr_t;

      return Expr<expr_t>(expr_t(expr1, expr2));
    }

    template <typename T1, typename T2>
    inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, ACosQuadOp > >
    acos_quad (const Expr<T1>& expr1, const Expr<T2>& expr2)
    {
      typedef BinaryExpr< Expr<T1>, Expr<T2>, ACosQuadOp > expr_t;

      return Expr<expr_t>(expr_t(expr1, expr2));
    }

    template <typename T1, typename T2>
    inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, ATanQuadOp > >
    atan_quad (const Expr<T1>& expr1, const Expr<T2>& expr2)
    {
      typedef BinaryExpr< Expr<T1>, Expr<T2>, ATanQuadOp > expr_t;

      return Expr<expr_t>(expr_t(expr1, expr2));
    }

    template <typename ExprT1, typename ExprT2>
    struct PowExprType {
      typedef UnaryExpr< ExprT1, LogOp > T3;
      typedef BinaryExpr< ExprT2, Expr<T3>, MultiplicationOp > T4;
      typedef UnaryExpr< Expr<T4>, ExpOp > T5;

      typedef Expr<T5> expr_type;
    };

     template <typename ExprT2>
     struct PowExprType< typename ExprT2::value_type, ExprT2 > {
       typedef typename ExprT2::value_type T1;
      typedef BinaryExpr< ExprT2, ConstExpr<T1>, MultiplicationOp > T4;
      typedef UnaryExpr< Expr<T4>, ExpOp > T5;

      typedef Expr<T5> expr_type;
    };

    template <typename ExprT1>
    struct PowExprType<ExprT1,typename ExprT1::value_type> {
      typedef typename ExprT1::value_type T2;
      typedef UnaryExpr< ExprT1, LogOp > T3;
      typedef BinaryExpr< ConstExpr<T2>, Expr<T3>, MultiplicationOp > T4;
      typedef UnaryExpr< Expr<T4>, ExpOp > T5;

      typedef Expr<T5> expr_type;
    };

    template <typename T1, typename T2>
    inline typename PowExprType< Expr<T1>, Expr<T2> >::expr_type
    pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
    {
      // pow(x,y) = exp(y*log(x))
      return exp(expr2*log(expr1));
    }

    template <typename T>
    inline typename PowExprType< typename Expr<T>::value_type, Expr<T> >::expr_type
    pow (const typename Expr<T>::value_type& c, const Expr<T>& expr)
    {
      // pow(x,y) = exp(y*log(x))
      return exp(expr*std::log(c));
    }

    template <typename T>
    inline typename PowExprType< Expr<T>, typename Expr<T>::value_type >::expr_type
    pow (const Expr<T>& expr, const typename Expr<T>::value_type& c)
    {
      // pow(x,y) = exp(y*log(x))
      return exp(c*log(expr));
    }

    template <typename T>
    struct Log10ExprType {
      typedef UnaryExpr< Expr<T>, LogOp > T1;
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;
      typedef BinaryExpr< Expr<T1>, ConstT, DivisionOp > T2;
      typedef Expr<T2> expr_type;
    };

    template <typename T>
    inline typename Log10ExprType<T>::expr_type
    log10 (const Expr<T>& expr)
    {
      // log10(x) = log(x)/log(10)
      return log(expr)/std::log(typename Expr<T>::value_type(10));
    }

    template <typename T>
    struct TanExprType {
      typedef UnaryExpr< Expr<T>, SinOp > T1;
      typedef UnaryExpr< Expr<T>, CosOp > T2;
      typedef BinaryExpr< Expr<T1>, Expr<T2>, DivisionOp > T3;
      typedef Expr<T3> expr_type;
    };

    template <typename T>
    inline typename TanExprType<T>::expr_type
    tan (const Expr<T>& expr)
    {
      // tan(x) = sin(x)/cos(x)
      return sin(expr)/cos(expr);
    }

    template <typename T>
    struct ASinExprType {
      typedef BinaryExpr< Expr<T>, Expr<T>, MultiplicationOp > T1;
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;
      typedef BinaryExpr< ConstT, Expr<T1>, SubtractionOp > T2;
      typedef UnaryExpr< Expr<T2>, SqrtOp > T3;
      typedef BinaryExpr< ConstT, Expr<T3>, DivisionOp > T4;
      typedef BinaryExpr< Expr<T>, Expr<T4>, ASinQuadOp > T5;
      typedef Expr<T5> expr_type;
    };

    template <typename T>
    inline typename ASinExprType<T>::expr_type
    asin (const Expr<T>& expr)
    {
      typedef typename Expr<T>::value_type value_type;

      // asin(x) = integral of 1/sqrt(1-x*x)
      return asin_quad(expr, value_type(1)/sqrt(value_type(1)-expr*expr));
    }

    template <typename T>
    struct ACosExprType {
      typedef BinaryExpr< Expr<T>, Expr<T>, MultiplicationOp > T1;
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;
      typedef BinaryExpr< ConstT, Expr<T1>, SubtractionOp > T2;
      typedef UnaryExpr< Expr<T2>, SqrtOp > T3;
      typedef BinaryExpr< ConstT, Expr<T3>, DivisionOp > T4;
      typedef BinaryExpr< Expr<T>, Expr<T4>, ACosQuadOp > T5;
      typedef Expr<T5> expr_type;
    };

    template <typename T>
    inline typename ACosExprType<T>::expr_type
    acos (const Expr<T>& expr)
    {
      typedef typename Expr<T>::value_type value_type;

      // acos(x) = integral of -1/sqrt(1-x*x)
      return acos_quad(expr, value_type(-1)/sqrt(value_type(1)-expr*expr));
    }

    template <typename T>
    struct ATanExprType {
      typedef BinaryExpr< Expr<T>, Expr<T>, MultiplicationOp > T1;
      typedef ConstExpr<typename Expr<T>::value_type> ConstT;
      typedef BinaryExpr< ConstT, Expr<T1>, AdditionOp > T2;
      typedef BinaryExpr< ConstT, Expr<T2>, DivisionOp > T3;
      typedef BinaryExpr< Expr<T>, Expr<T3>, ATanQuadOp > T4;
      typedef Expr<T4> expr_type;
    };

    template <typename T>
    inline typename ATanExprType<T>::expr_type
    atan (const Expr<T>& expr)
    {
      typedef typename Expr<T>::value_type value_type;

      // atan(x) = integral of 1/(1+x*x)
      return atan_quad(expr, value_type(1)/(value_type(1)+expr*expr));
    }

    template <typename T>
    struct TanhExprType {
      typedef UnaryExpr< Expr<T>, SinhOp > T1;
      typedef UnaryExpr< Expr<T>, CoshOp > T2;
      typedef BinaryExpr< Expr<T1>, Expr<T2>, DivisionOp > T3;
      typedef Expr<T3> expr_type;
    };

    template <typename T>
    inline typename TanhExprType<T>::expr_type
    tanh (const Expr<T>& expr)
    {
      // tanh(x) = sinh(x)/cosh(x)
      return sinh(expr)/cosh(expr);
    }

  } // namespace Tay

} // namespace Sacado

namespace std {
  using Sacado::Tay::pow;
  using Sacado::Tay::log10;
  using Sacado::Tay::tan;
  using Sacado::Tay::asin;
  using Sacado::Tay::acos;
  using Sacado::Tay::atan;
  using Sacado::Tay::tanh;
}

//-------------------------- Relational Operators -----------------------

#define TAYLOR_RELOP_MACRO(OP)                                          \
namespace Sacado {                                                      \
  namespace Tay {                                                       \
    template <typename ExprT1, typename ExprT2>                         \
    inline bool                                                         \
    operator OP (const Expr<ExprT1>& expr1,                             \
                 const Expr<ExprT2>& expr2)                             \
    {                                                                   \
      return expr1.fastAccessCoeff(0) OP expr2.fastAccessCoeff(0);      \
    }                                                                   \
                                                                        \
    template <typename ExprT2>                                          \
    inline bool                                                         \
    operator OP (const typename Expr<ExprT2>::value_type& a,            \
                 const Expr<ExprT2>& expr2)                             \
    {                                                                   \
      return a OP expr2.fastAccessCoeff(0);                             \
    }                                                                   \
                                                                        \
    template <typename ExprT1>                                          \
    inline bool                                                         \
    operator OP (const Expr<ExprT1>& expr1,                             \
                 const typename Expr<ExprT1>::value_type& b)            \
    {                                                                   \
      return expr1.fastAccessCoeff(0) OP b;                             \
    }                                                                   \
  }                                                                     \
}

TAYLOR_RELOP_MACRO(==)
TAYLOR_RELOP_MACRO(!=)
TAYLOR_RELOP_MACRO(<)
TAYLOR_RELOP_MACRO(>)
TAYLOR_RELOP_MACRO(<=)
TAYLOR_RELOP_MACRO(>=)
TAYLOR_RELOP_MACRO(<<=)
TAYLOR_RELOP_MACRO(>>=)
TAYLOR_RELOP_MACRO(&)
TAYLOR_RELOP_MACRO(|)

#undef TAYLOR_RELOP_MACRO

namespace Sacado {

  namespace Tay {

    template <typename ExprT>
    inline bool operator ! (const Expr<ExprT>& expr)
    {
      return ! expr.fastAccessCoeff(0);
    }

  } // namespace Tay

} // namespace Sacado

//-------------------------- Boolean Operators -----------------------
namespace Sacado {

  namespace Tay {

    template <typename ExprT>
    bool toBool2(const Expr<ExprT>& x) {
      bool is_zero = true;
      for (int i=0; i<=x.degree(); i++)
        is_zero = is_zero && (x.coeff(i) == 0.0);
      return !is_zero;
    }

  } // namespace Tay

} // namespace Sacado

#define TAY_BOOL_MACRO(OP)                                              \
namespace Sacado {                                                      \
  namespace Tay {                                                       \
    template <typename ExprT1, typename ExprT2>                         \
    inline bool                                                         \
    operator OP (const Expr<ExprT1>& expr1,                             \
                 const Expr<ExprT2>& expr2)                             \
    {                                                                   \
      return toBool2(expr1) OP toBool2(expr2);                          \
    }                                                                   \
                                                                        \
    template <typename ExprT2>                                          \
    inline bool                                                         \
    operator OP (const typename Expr<ExprT2>::value_type& a,            \
                 const Expr<ExprT2>& expr2)                             \
    {                                                                   \
      return a OP toBool2(expr2);                                       \
    }                                                                   \
                                                                        \
    template <typename ExprT1>                                          \
    inline bool                                                         \
    operator OP (const Expr<ExprT1>& expr1,                             \
                 const typename Expr<ExprT1>::value_type& b)            \
    {                                                                   \
      return toBool2(expr1) OP b;                                       \
    }                                                                   \
  }                                                                     \
}

TAY_BOOL_MACRO(&&)
TAY_BOOL_MACRO(||)

#undef TAY_BOOL_MACRO

//-------------------------- I/O Operators -----------------------

namespace Sacado {

  namespace Tay {

    template <typename ExprT>
    std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
      os.setf(std::ios::fixed, std::ios::floatfield);
      os.width(12);
      os << "[";

      for (int i=0; i<=x.degree(); i++) {
        os.width(12);
        os << x.coeff(i);
      }

      os << "]";
      return os;
    }

    template <typename T>
    CacheTaylor<T> diff(const CacheTaylor<T>& x, int n = 1) {
      const int d = x.degree();
      if (n <= 0)
        return x;
      else if (n > d) {
        Taylor<T> y(0);
        return y;
      }
      CacheTaylor<T> y(d-n);
      int c = 1;
      for (int i=1; i<=n; ++i)
        c *= i;
      for (int i=n; i<=d; ++i) {
        y.fastAccessCoeff(i-n) = x.fastAccessCoeff(i) * T(c);
        c = (c / (i-n+1)) * (i+1);
      }
      return y;
    }

  } // namespace Tay

} // namespace Sacado


#endif // SACADO_TAY_CACHETAYLOROPS_HPP