Sacado_Tay_Taylor.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                           Sacado Package
//                 Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// 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_TAYLOR_HPP
#define SACADO_TAY_TAYLOR_HPP

#include "Sacado_ConfigDefs.h"
#include "Sacado_Base.hpp"
#include "Sacado_Handle.hpp"
#include <cmath>
#include <algorithm>    // for std::min and std::max
#include <ostream>      // for std::ostream
#include "Sacado_dummy_arg.hpp"

namespace Sacado {

  namespace Tay {


    template <typename T>
    class Taylor : public Base< Taylor<T> > {
    public:

      template <typename U>
      struct apply {
        typedef Taylor<U> type;
      };

      typedef T value_type;

      typedef typename ScalarType<value_type>::type scalar_type;

      Taylor();


      Taylor(const T& x);


      Taylor(const typename dummy<value_type,scalar_type>::type& x);


      Taylor(int d, const T & x);


      explicit Taylor(int d);

      Taylor(const Taylor& x);

      ~Taylor();


      void resize(int d, bool keep_coeffs);


      void reserve(int d);


      void copyForWrite() { th.makeOwnCopy(); }

      bool isEqualTo(const Taylor& x) const {
        typedef IsEqual<value_type> IE;
        if (x.degree() != this->degree()) return false;
        bool eq = true;
        for (int i=0; i<=this->degree(); i++)
          eq = eq && IE::eval(x.coeff(i), this->coeff(i));
        return eq;
      }


      Taylor<T>& operator=(const T& val);


      Taylor<T>& operator=(const typename dummy<value_type,scalar_type>::type& val);

      Taylor<T>& operator=(const Taylor<T>& x);



      const T& val() const { return th->coeff_[0];}

      T& val() { return th->coeff_[0];}



      int degree() const { return th->deg_;}

      bool hasFastAccess(int d) const { return th->deg_>=d;}

      const T* coeff() const { return th->coeff_;}

      T* coeff() { return th->coeff_;}

      T coeff(int i) const {
        T tmp= i<=th->deg_ ? th->coeff_[i]:T(0.); return tmp;}

      T& fastAccessCoeff(int i) { return th->coeff_[i];}

      const T& fastAccessCoeff(int i) const { return th->coeff_[i];}



      Taylor<T> operator + () const;

      Taylor<T> operator - () const;

      Taylor<T>& operator += (const T& x);

      Taylor<T>& operator -= (const T& x);

      Taylor<T>& operator *= (const T& x);

      Taylor<T>& operator /= (const T& x);

      Taylor<T>& operator += (const Taylor<T>& x);

      Taylor<T>& operator -= (const Taylor<T>& x);

      Taylor<T>& operator *= (const Taylor<T>& x);

      Taylor<T>& operator /= (const Taylor<T>& x);


    protected:

      int length() const { return th->len_; }

      void resizeCoeffs(int len);

    protected:

      struct TaylorData {

        T* coeff_;

        int deg_;

        int len_;

        TaylorData();

        TaylorData(const T& x);

        TaylorData(int d, const T & x);

        TaylorData(int d);

        TaylorData(int d, int l);

        TaylorData(const TaylorData& x);

        ~TaylorData();

        TaylorData& operator=(const TaylorData& x);

      };

      Sacado::Handle<TaylorData> th;

    }; // class Taylor

    template <typename T>
    Taylor<T> diff(const Taylor<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;
      }
      Taylor<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;
    }

    // Operations
    template <typename T> Taylor<T> operator+(const Base< Taylor<T> >& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator+(const typename Taylor<T>::value_type& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator+(const Base< Taylor<T> >& a,
                                              const typename Taylor<T>::value_type& b);
    template <typename T> Taylor<T> operator-(const Base< Taylor<T> >& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator-(const typename Taylor<T>::value_type& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator-(const Base< Taylor<T> >& a,
                                              const typename Taylor<T>::value_type& b);
    template <typename T> Taylor<T> operator*(const Base< Taylor<T> >& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator*(const typename Taylor<T>::value_type& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator*(const Base< Taylor<T> >& a,
                                              const typename Taylor<T>::value_type& b);
    template <typename T> Taylor<T> operator/(const Base< Taylor<T> >& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator/(const typename Taylor<T>::value_type& a,
                                              const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> operator/(const Base< Taylor<T> >& a,
                                              const typename Taylor<T>::value_type& b);
    template <typename T> Taylor<T> exp(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> log(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> log10(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> sqrt(const Base< Taylor<T> >& a);
#ifdef HAVE_SACADO_CXX11
    template <typename T> Taylor<T> cbrt(const Base< Taylor<T> >& a);
#endif
    template <typename T> Taylor<T> pow(const Base< Taylor<T> >& a,
                                        const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> pow(const typename Taylor<T>::value_type& a,
                                        const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> pow(const Base< Taylor<T> >& a,
                                        const typename Taylor<T>::value_type& b);
    template <typename T> void sincos(const Base< Taylor<T> >& a,
                                      Taylor<T>& s, Taylor<T>& c);
    template <typename T> Taylor<T> cos(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> sin(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> tan(const Base< Taylor<T> >& a);
    template <typename T> void sinhcosh(const Base< Taylor<T> >& a,
                                        Taylor<T>& s, Taylor<T>& c);
    template <typename T> Taylor<T> cosh(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> sinh(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> tanh(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> quad(const typename Taylor<T>::value_type& c0,
                                         const Base< Taylor<T> >& a,
                                         const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> acos(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> asin(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> atan(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> atan2(const Base< Taylor<T> >& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> atan2(const typename Taylor<T>::value_type& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> atan2(const Base< Taylor<T> >& a,
                                          const typename Taylor<T>::value_type& b);
    template <typename T> Taylor<T> acosh(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> asinh(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> atanh(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> abs(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> fabs(const Base< Taylor<T> >& a);
    template <typename T> Taylor<T> max(const Base< Taylor<T> >& a,
                                        const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> max(const typename Taylor<T>::value_type& a,
                                        const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> max(const Base< Taylor<T> >& a,
                                        const typename Taylor<T>::value_type& b);
    template <typename T> Taylor<T> min(const Base< Taylor<T> >& a,
                                        const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> min(const typename Taylor<T>::value_type& a,
                                        const Base< Taylor<T> >& b);
    template <typename T> Taylor<T> min(const Base< Taylor<T> >& a,
                                        const typename Taylor<T>::value_type& b);
    template <typename T> bool operator==(const Base< Taylor<T> >& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator==(const typename Taylor<T>::value_type& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator==(const Base< Taylor<T> >& a,
                                          const typename Taylor<T>::value_type& b);
    template <typename T> bool operator!=(const Base< Taylor<T> >& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator!=(const typename Taylor<T>::value_type& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator!=(const Base< Taylor<T> >& a,
                                          const typename Taylor<T>::value_type& b);
    template <typename T> bool operator<=(const Base< Taylor<T> >& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator<=(const typename Taylor<T>::value_type& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator<=(const Base< Taylor<T> >& a,
                                          const typename Taylor<T>::value_type& b);
    template <typename T> bool operator>=(const Base< Taylor<T> >& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator>=(const typename Taylor<T>::value_type& a,
                                          const Base< Taylor<T> >& b);
    template <typename T> bool operator>=(const Base< Taylor<T> >& a,
                                          const typename Taylor<T>::value_type& b);
    template <typename T> bool operator<(const Base< Taylor<T> >& a,
                                         const Base< Taylor<T> >& b);
    template <typename T> bool operator<(const typename Taylor<T>::value_type& a,
                                         const Base< Taylor<T> >& b);
    template <typename T> bool operator<(const Base< Taylor<T> >& a,
                                         const typename Taylor<T>::value_type& b);
    template <typename T> bool operator>(const Base< Taylor<T> >& a,
                                         const Base< Taylor<T> >& b);
    template <typename T> bool operator>(const typename Taylor<T>::value_type& a,
                                         const Base< Taylor<T> >& b);
    template <typename T> bool operator>(const Base< Taylor<T> >& a,
                                         const typename Taylor<T>::value_type& b);
    template <typename T> std::ostream& operator << (std::ostream& os,
                                                     const Base< Taylor<T> >& a);

  } // namespace Tay

} // namespace Sacado

#include "Sacado_Tay_TaylorTraits.hpp"
#include "Sacado_Tay_TaylorImp.hpp"

#endif // SACADO_TAY_TAYLOR_HPP