JMathlib.hh

Go to the documentation of this file.

#ifndef __JMATH__JMATHLIB__
#define __JMATH__JMATHLIB__
 
#include <cstddef>
#include <cmath>
#include <array>
#include <utility>
#include <memory>
 
 
/**
 * \file
 * Functional algebra.
 *
 * \author mdejong
 */
 
namespace JMATH {}
namespace JPP { using namespace JMATH; }
 
namespace JMATH {
 
 
  /**
   * Type definition for list of parameters.
   *
   * The template argument refers to the function that could be fitted to some data.\n
   * The parameters of the function should be data members of type <tt>double</tt>.
   */
  template<class JF1_t, size_t N>
  using parameter_list = std::array<double JF1_t::*, N>;
 
 
  /**
   * Auxiliary method to join parameter and list of parameters.
   *
   * \param  parameters       parameter list
   * \param  ls               indices of parameter list
   * \param  parameter        parameter
   * \return                  parameter list
   */
  template<typename U, typename V, typename W, size_t ...i>
  constexpr parameter_list<U, sizeof...(i) + 1> join(const parameter_list<V, sizeof...(i)>& parameters,
                                                     const std::index_sequence<i...>& ls,
                                                     double W::* parameter)
  {
        return {parameters[i]..., parameter};
  };
 
 
  /**
   * Join parameter and list of parameters.
   *
   * \param  parameters       parameter list
   * \param  parameter        parameter
   * \return                  parameter list
   */
  template<typename U, typename V, size_t N, typename W>
  constexpr parameter_list<U, N + 1> join(const parameter_list<V, N>& parameters,
                                          double W::* parameter)
  {
        return join<U>(parameters, std::make_index_sequence<N>{}, parameter);
  }
 
 
  /**
   * Auxiliary method to join two parameter lists.
   *
   * \param  v                first  parameter list
   * \param  lv               indices of first  parameter list
   * \param  w                second parameter list
   * \param  lw               indices of second parameter list
   * \return                  parameter list
   */
  template<typename U, typename V, typename W, size_t ...iv, size_t ...iw>
  constexpr parameter_list<U, sizeof...(iv) + sizeof...(iw)> join(const parameter_list<V, sizeof...(iv)>& v,
                                                                  const std::index_sequence<iv...>& lv,
                                                                  const parameter_list<W, sizeof...(iw)>& w,
                                                                  const std::index_sequence<iw...>& lw)
  {
        return {v[iv]..., w[iw]... };
  };
 
 
  /**
   * Join two parameter lists.
   *
   * \param  v                first  parameter list
   * \param  w                second parameter list
   * \return                  parameter list
   */
  template<typename U, typename V, size_t NV, typename W, size_t NW>
  constexpr parameter_list<U, NV + NW> join(const parameter_list<V, NV>& v,
                                            const parameter_list<W, NW>& w)
  {
        return join<U>(v, std::make_index_sequence<NV>{},
                   w, std::make_index_sequence<NW>{});
  }
 
 
  /**
   * Auxiliary method to convert parameter list.
   *
   * \param  parameters       parameter list
   * \param  ls               indices of parameter list
   * \return                  parameter list
   */
  template<typename U, typename V, size_t ...i>
  constexpr parameter_list<U, sizeof...(i)> convert(const parameter_list<V, sizeof...(i)>& parameters,
                                                    const std::index_sequence<i...>& ls)
  {
        return {parameters[i]...};
  };
 
 
  /**
   * Convert parameter list.
   *
   * \param  parameters       parameter list
   * \return                  parameter list
   */
  template<typename U, typename V, size_t N>
  constexpr parameter_list<U, N> convert(const parameter_list<V, N>& parameters)
  {
        return convert<U>(parameters, std::make_index_sequence<N>{});
  }
 
 
  /**
   * Get number of parameters.
   *
   * \return               number of parameters
   */
  template<class JF1_t>
  constexpr size_t getNumberOfParameters()
  {
    return JF1_t::parameters.size();
  }
 
 
  /**
   * Set values of all parameters.
   *
   * \param  f1            pointer to function
   * \param  values        pointer to list of values
   */
  template<class JF1_t>
  inline void setParameters(JF1_t* f1, const double* values)
  {
    if (f1 != NULL && values != NULL) {
      for (size_t i = 0; i != getNumberOfParameters<JF1_t>(); ++i) {
        (*f1)[i] = values[i];
      }
    }
  }
 
 
  /**
   * Auxiliary base class for mathematical operations on parameters of function.
   */
  template<class JF1_t>
  struct JCalculus 
  {
    /**
     * Negate function.
     *
     * \return               this function
     */
    JF1_t& negate()
    {
      for (const auto& i : JF1_t::parameters) {
        static_cast<JF1_t&>(*this).*i = -(static_cast<JF1_t&>(*this).*i);
      }
 
      return static_cast<JF1_t&>(*this);
    }
 
 
    /**
     * Add function.
     *
     * \param  f1            function
     * \return               this function
     */
    JF1_t& add(const JF1_t& f1)
    {
      for (const auto& i : JF1_t::parameters) {
        static_cast<JF1_t&>(*this).*i += f1.*i;
      }
 
      return static_cast<JF1_t&>(*this);
    }
 
 
    /**
     * Subtract function.
     *
     * \param  f1            function
     * \return               this function
     */
    JF1_t& sub(const JF1_t& f1)
    {
      for (const auto& i : JF1_t::parameters) {
        static_cast<JF1_t&>(*this).*i -= f1.*i;
      }
 
      return static_cast<JF1_t&>(*this);
    }
 
 
    /**
     * Scale function.
     *
     * \param  factor        factor
     * \return               this function
     */
    JF1_t& mul(const double factor)
    {
      for (const auto& i : JF1_t::parameters) {
        static_cast<JF1_t&>(*this).*i *= factor;
      }
 
      return static_cast<JF1_t&>(*this);
    }
 
 
    /**
     * Scale function.
     *
     * \param  factor        factor
     * \return               this function
     */
    JF1_t& div(const double factor)
    {
      for (const auto& i : JF1_t::parameters) {
        static_cast<JF1_t&>(*this).*i /= factor;
      }
 
      return static_cast<JF1_t&>(*this);
    }
 
 
    /**
     * Add function.
     *
     * \param  function       this function
     * \param  value          value
     * \return                this function
     */
    friend JF1_t& operator+=(JF1_t& function, const JF1_t& value)
    {
      return function.add(value);
    }
 
 
    /**
     * Subtract function.
     *
     * \param  function       this function
     * \param  value          value
     * \return                this function
     */
    friend JF1_t& operator-=(JF1_t& function, const JF1_t& value)
    {
      return function.sub(value);
    }
 
 
    /**
     * Scale function.
     *
     * \param  function       this function
     * \param  factor         factor
     * \return                this function
     */
    friend JF1_t& operator*=(JF1_t& function, const double factor)
    {
      return function.mul(factor);
    }
 
 
    /**
     * Scale function.
     *
     * \param  function       this function
     * \param  factor         factor
     * \return                this function
     */
    friend JF1_t& operator/=(JF1_t& function, const double factor)
    {
      return function.div(factor);
    }
  };
 
 
  template<class JF1_t>
  struct JNegate;            // forward declaration for negation of function.
 
  template<class JF1_t, class JF2_t = double>
  struct JAdd;               // forward declaration for addition of fuction.
 
  template<class JF1_t, class JF2_t = double>
  struct JSub;               // forward declaration for subtraction of fuction.
 
  template<class JF1_t, class JF2_t = double>
  struct JMul;               // forward declaration for multiplication of fuction.
 
  template<class JF1_t, class JF2_t = double>
  struct JDiv;               // forward declaration for division of fuction.
 
  template<class JF1_t = void>
  struct JFn;                // forward declaration for fixed power of function.
 
  template<class JF1_t = void>
  struct JPf;                // forward declaration for pointer to function.
 
 
  /**
   * Auxiliary base class for mathematical operations on functions.
   */
  template<class JF1_t>
  struct JMathlib {
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double operator()(const Args& ...args) const
    {
      return static_cast<const JF1_t&>(*this).getValue(args...);
    }
 
 
    /**
     * Affirm operator.
     *
     * \param  function      this function
     * \return               result function
     */
    friend const JF1_t& operator+(const JF1_t& function)
    {
      return function;
    }
 
 
    /**
     * Negate operator.
     *
     * \param  function      this function
     * \return               result function
     */
    friend JNegate<JF1_t> operator-(const JF1_t& function)
    {
      return JNegate<JF1_t>(function);
    }
 
 
    /**
     * Addition of constant value.
     *
     * \param  f1            function
     * \param  value         value
     * \return               result function
     */
    friend JAdd<JF1_t> operator+(const JF1_t& f1, const double value)
    {
      return JAdd<JF1_t>(f1, value);
    }
 
 
    /**
     * Addition of constant value.
     *
     * \param  value         value
     * \param  f1            function
     * \return               result function
     */
    friend JAdd<JF1_t> operator+(const double value, const JF1_t& f1)
    {
      return JAdd<JF1_t>(f1, value);
    }
 
 
    /**
     * Subtraction of constant value.
     *
     * \param  f1            function
     * \param  value         value
     * \return               result function
     */
    friend JSub<JF1_t> operator-(const JF1_t& f1, const double value)
    {
      return JSub<JF1_t>(f1, value);
    }
 
 
    /**
     * Subtraction of constant value.
     *
     * \param  f1            function
     * \param  value         value
     * \return               result function
     */
    friend JAdd< JNegate<JF1_t> > operator-(const double value, const JF1_t& f1)
    {
      return JAdd< JNegate<JF1_t> >(JNegate<JF1_t>(f1), value);
    }
 
 
    /**
     * Multiplication of constant value.
     *
     * \param  f1            function
     * \param  value         value
     * \return               result function
     */
    friend JMul<JF1_t> operator*(const JF1_t& f1, const double value)
    {
      return JMul<JF1_t>(f1, value);
    }
 
 
    /**
     * Multiplication of constant value.
     *
     * \param  value         value
     * \param  f1            function
     * \return               result function
     */
    friend JMul<JF1_t> operator*(const double value, const JF1_t& f1)
    {
      return JMul<JF1_t>(f1, value);
    }
 
 
    /**
     * Division of constant value.
     *
     * \param  f1            function
     * \param  value         value
     * \return               result function
     */
    friend JDiv<JF1_t> operator/(const JF1_t& f1, const double value)
    {
      return JDiv<JF1_t>(f1, value);
    }
 
 
    /**
     * Addition of two functions.
     *
     * \param  f1            first  function
     * \param  f2            second function
     * \return               result function
     */
    template<class JF2_t>
    friend JAdd<JF1_t, JF2_t> operator+(const JF1_t& f1, const JF2_t& f2)
    {
      return JAdd<JF1_t, JF2_t>(f1, f2);
    }
 
 
    /**
     * Subtraction of two functions.
     *
     * \param  f1            first  function
     * \param  f2            second function
     * \return               result function
     */
    template<class JF2_t>
    friend JSub<JF1_t, JF2_t> operator-(const JF1_t& f1, const JF2_t& f2)
    {
      return JSub<JF1_t, JF2_t>(f1, f2);
    }
 
 
    /**
     * Multiplication of two functions.
     *
     * \param  f1            first  function
     * \param  f2            second function
     * \return               result function
     */
    template<class JF2_t>
    friend JMul<JF1_t, JF2_t> operator*(const JF1_t& f1, const JF2_t& f2)
    {
      return JMul<JF1_t, JF2_t>(f1, f2);
    }
 
 
    /**
     * Division of two functions.
     *
     * \param  f1            first  function
     * \param  f2            second function
     * \return               result function
     */
    template<class JF2_t>
    friend JDiv<JF1_t, JF2_t> operator/(const JF1_t& f1, const JF2_t& f2)
    {
      return JDiv<JF1_t, JF2_t>(f1, f2);
    }
 
 
    /**
     * Power-of operator.
     *
     * \param  f1            function
     * \param  N             power
     * \return               result function
     */
    friend JFn<JF1_t> operator^(const JF1_t& f1, int N)
    {
      return JFn<JF1_t>(f1, N);
    }
 
 
    /**
     * Address-of operator.
     *
     * \param  f1            function
     * \return               result function
     */
    friend JPf<JF1_t> operator&(const JF1_t& f1)
    {
      return JPf<JF1_t>(std::addressof(f1));
    }
 
 
    /**
     * Get value of parameter at given index.
     *
     * \param  i             index
     * \return               value
     */
    double operator[](const size_t i) const
    {
      return static_cast<const JF1_t&>(*this).*JF1_t::parameters[i];
    }
 
 
    /**
     * Get value of parameter at given index.
     *
     * \param  i             index
     * \return               value
     */
    double& operator[](const size_t i)
    {
      return static_cast<JF1_t&>(*this).*JF1_t::parameters[i];
    }
  };
 
 
  /**
   * Negate of function.
   */
  template<class JF1_t>
  struct JNegate :
    public JMathlib< JNegate<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JNegate<JF1_t> >::operator();
    using JMathlib< JNegate<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JNegate()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JNegate(const JF1_t& f1) :
      JF1_t(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return -static_cast<const JF1_t&>(*this).getValue(args...);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return -static_cast<const JF1_t&>(*this).getDerivative(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JNegate getGradient(const Args& ...args) const
    {
      JNegate gradient;
 
      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF1_t&>(gradient).negate();
 
      return gradient;
    }
  };
 
 
  /**
   * Addition of constant value.
   */
  template<class JF1_t>
  struct JAdd<JF1_t, double> :
    public JMathlib< JAdd<JF1_t, double> >,
    public JF1_t
  {
    using JMathlib< JAdd<JF1_t, double> >::operator();
    using JMathlib< JAdd<JF1_t, double> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JAdd()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     * \param  value         value
     */
    JAdd(const JF1_t& f1, const double value) :
      JF1_t(f1),
      value(value)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return static_cast<const JF1_t&>(*this).getValue(args...) + value;
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return static_cast<const JF1_t&>(*this).getDerivative(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JAdd getGradient(const Args& ...args) const
    {
      JAdd gradient;
 
      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(args...);
 
      return gradient;
    }
 
  private:
    double value;
  };
 
 
  /**
   * Subtraction of constant value.
   */
  template<class JF1_t>
  struct JSub<JF1_t, double> :
    public JMathlib< JSub<JF1_t, double> >,
    public JF1_t
  {
    using JMathlib< JSub<JF1_t, double> >::operator();
    using JMathlib< JSub<JF1_t, double> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JSub()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     * \param  value         value
     */
    JSub(const JF1_t& f1, const double value) :
      JF1_t(f1),
      value(value)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return static_cast<const JF1_t&>(*this).getValue(args...) - value;
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return static_cast<const JF1_t&>(*this).getDerivative(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JSub getGradient(const Args& ...args) const
    {
      JSub gradient;
 
      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(args...);
 
      return gradient;
    }
 
  private:
    double value;
  };
 
 
  /**
   * Multiplication of constant value.
   */
  template<class JF1_t>
  struct JMul<JF1_t, double> :
    public JMathlib< JMul<JF1_t, double> >,
    public JF1_t
  {
    using JMathlib< JMul<JF1_t, double> >::operator();
    using JMathlib< JMul<JF1_t, double> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JMul()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     * \param  value         value
     */
    JMul(const JF1_t& f1, const double value) :
      JF1_t(f1),
      value(value)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return static_cast<const JF1_t&>(*this).getValue(args...) * value;
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return static_cast<const JF1_t&>(*this).getDerivative(x) * value;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JMul getGradient(const Args& ...args) const
    {
      JMul gradient;
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF1_t&>(gradient) *= value;
 
      return gradient;
    }
 
  private:
    double value;
  };
 
 
  /**
   * Division of constant value.
   */
  template<class JF1_t>
  struct JDiv<JF1_t, double> :
    public JMathlib< JDiv<JF1_t, double> >,
    public JF1_t
  {
    using JMathlib< JDiv<JF1_t, double> >::operator();
    using JMathlib< JDiv<JF1_t, double> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JDiv()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     * \param  value         value
     */
    JDiv(const JF1_t& f1, const double value) :
      JF1_t(f1),
      value(value)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return static_cast<const JF1_t&>(*this).getValue(args...) / value;
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return static_cast<const JF1_t&>(*this).getDerivative(x) / value;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JDiv getGradient(const Args& ...args) const
    {
      JDiv gradient;
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF1_t&>(gradient) /= value;
 
      return gradient;
    }
 
  private:
    double value;
  };
 
 
  /**
   * Auxiliary data structure for pair of functions.
   */
  template<class JF1_t, class JF2_t,
           bool v1 = JF1_t::parameters.size() == 0 && !std::is_base_of<JPf<>, JF1_t>::value && !std::is_base_of<JFn<>, JF1_t>::value,   // true -> virtual base class
           bool v2 = JF2_t::parameters.size() == 0 && !std::is_base_of<JPf<>, JF2_t>::value && !std::is_base_of<JFn<>, JF2_t>::value>   // true -> virtual base class
  struct JPair :
    public JF1_t,
    public JF2_t
  {
    /**
     * Default constructor.
     */
    JPair()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JPair(const JF1_t& f1, const JF2_t& f2) :
      JF1_t(f1),
      JF2_t(f2)
    {}
 
    static constexpr parameter_list<JPair, JF1_t::parameters.size() + JF2_t::parameters.size()> parameters = join<JPair>(JF1_t::parameters, JF2_t::parameters);
  };
 
 
  /**
   * Template specialisation for pair of functions.
   */
  template<class JF1_t, class JF2_t>
  struct JPair<JF1_t, JF2_t, false, true> :
    public virtual JF2_t,
    public JF1_t
  {
    /**
     * Default constructor.
     */
    JPair()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JPair(const JF1_t& f1, const JF2_t& f2) :
      JF1_t(f1)
    {}
 
    static constexpr parameter_list<JPair, JF1_t::parameters.size()> parameters = convert<JPair>(JF1_t::parameters);
  };
 
 
  /**
   * Template specialisation for pair of functions.
   */
  template<class JF1_t, class JF2_t>
  struct JPair<JF1_t, JF2_t, true, false> :
    public virtual JF1_t,
    public JF2_t
  {
    /**
     * Default constructor.
     */
    JPair()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JPair(const JF1_t& f1, const JF2_t& f2) :
      JF2_t(f2)
    {}
 
    static constexpr parameter_list<JPair, JF2_t::parameters.size()> parameters = convert<JPair>(JF2_t::parameters);
  };
 
 
  /**
   * Template specialisation for pair of functions.
   */
  template<class JF1_t, class JF2_t>
  struct JPair<JF1_t, JF2_t, true, true> :
    public virtual JF1_t,
    public virtual JF2_t
  {
    /**
     * Default constructor.
     */
    JPair()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JPair(const JF1_t& f1, const JF2_t& f2)
    {}
 
    static constexpr parameter_list<JPair, 0> parameters{};
  };
 
 
  /**
   * Auxiliary data structure for pair of functions.
   */
  template<class JF1_t, class JF2_t>
  struct JPair_t :
    public JPair<JF1_t, JF2_t>,
    public JCalculus< JPair_t<JF1_t, JF2_t> >
  {
    using JCalculus< JPair_t<JF1_t, JF2_t> >::negate;
    using JCalculus< JPair_t<JF1_t, JF2_t> >::add;
    using JCalculus< JPair_t<JF1_t, JF2_t> >::sub;
    using JCalculus< JPair_t<JF1_t, JF2_t> >::mul;
    using JCalculus< JPair_t<JF1_t, JF2_t> >::div;
 
 
    /**
     * Default constructor.
     */
    JPair_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JPair_t(const JF1_t& f1, const JF2_t& f2) :
      JPair<JF1_t, JF2_t>(f1, f2)
    {}
  };
 
 
  /**
   * Addition of two functions.
   */
  template<class JF1_t, class JF2_t>
  struct JAdd :
    public JMathlib< JAdd<JF1_t, JF2_t> >,
    public JPair_t<JF1_t, JF2_t>
  {
    using JMathlib< JAdd<JF1_t, JF2_t> >::operator();
    using JMathlib< JAdd<JF1_t, JF2_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JAdd()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JAdd(const JF1_t& f1, const JF2_t& f2) :
      JPair_t<JF1_t, JF2_t>(f1, f2)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return (static_cast<const JF1_t&>(*this).getValue(args...)  + 
              static_cast<const JF2_t&>(*this).getValue(args...));
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return (static_cast<const JF1_t&>(*this).getDerivative(x)  + 
              static_cast<const JF2_t&>(*this).getDerivative(x));
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JAdd getGradient(const Args& ...args) const
    {
      JAdd gradient;
 
      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF2_t&>(gradient) = static_cast<const JF2_t&>(*this).getGradient(args...);
 
      return gradient;
    }
  };
 
 
  /**
   * Subtraction of two functions.
   */
  template<class JF1_t, class JF2_t>
  struct JSub :
    public JMathlib< JSub<JF1_t, JF2_t> >,
    public JPair_t<JF1_t, JF2_t>
  {
    using JMathlib< JSub<JF1_t, JF2_t> >::operator();
    using JMathlib< JSub<JF1_t, JF2_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JSub()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JSub(const JF1_t& f1, const JF2_t& f2) :
      JPair_t<JF1_t, JF2_t>(f1, f2)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return (static_cast<const JF1_t&>(*this).getValue(args...)  - 
              static_cast<const JF2_t&>(*this).getValue(args...));
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return (static_cast<const JF1_t&>(*this).getDerivative(x)  - 
              static_cast<const JF2_t&>(*this).getDerivative(x));
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JSub getGradient(const Args& ...args) const
    {
      JSub gradient;
 
      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF2_t&>(gradient) = static_cast<const JF2_t&>(*this).getGradient(args...);
      static_cast<JF2_t&>(gradient).negate();
 
      return gradient;
    }
  };
 
 
  /**
   * Multiplication of two functions.
   */
  template<class JF1_t, class JF2_t>
  struct JMul :
    public JMathlib< JMul<JF1_t, JF2_t> >,
    public JPair_t<JF1_t, JF2_t>
  {
    using JMathlib< JMul<JF1_t, JF2_t> >::operator();
    using JMathlib< JMul<JF1_t, JF2_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JMul()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JMul(const JF1_t& f1, const JF2_t& f2) :
      JPair_t<JF1_t, JF2_t>(f1, f2)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return (static_cast<const JF1_t&>(*this).getValue(args...)  *
              static_cast<const JF2_t&>(*this).getValue(args...));
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return (static_cast<const JF1_t&>(*this).getDerivative(x) * static_cast<const JF2_t&>(*this).getValue(x)  +
              static_cast<const JF1_t&>(*this).getValue(x) * static_cast<const JF2_t&>(*this).getDerivative(x));
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JMul getGradient(const Args& ...args) const
    {
      JMul gradient;
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF1_t&>(gradient) *= static_cast<const JF2_t&>(*this).getValue(args...);
      static_cast<JF2_t&>(gradient)  = static_cast<const JF2_t&>(*this).getGradient(args...);
      static_cast<JF2_t&>(gradient) *= static_cast<const JF1_t&>(*this).getValue(args...);
 
      return gradient;
    }
  };
 
 
  /**
   * Division of two functions.
   */
  template<class JF1_t, class JF2_t>
  struct JDiv :
    public JMathlib< JDiv<JF1_t, JF2_t> >,
    public JPair_t<JF1_t, JF2_t>
  {
    using JMathlib< JDiv<JF1_t, JF2_t> >::operator();
    using JMathlib< JDiv<JF1_t, JF2_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JDiv()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  f2            second function
     */
    JDiv(const JF1_t& f1, const JF2_t& f2) :
      JPair_t<JF1_t, JF2_t>(f1, f2)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return (static_cast<const JF1_t&>(*this).getValue(args...)  /
              static_cast<const JF2_t&>(*this).getValue(args...));
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double v = static_cast<const JF1_t&>(*this).getValue(x);
      const double w = static_cast<const JF2_t&>(*this).getValue(x);
 
      return (static_cast<const JF1_t&>(*this).getDerivative(x) * w  - 
              v * static_cast<const JF2_t&>(*this).getDerivative(x)) / (w*w);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JDiv getGradient(const Args& ...args) const
    {
      JDiv gradient;
 
      const double v = static_cast<const JF1_t&>(*this).getValue(args...);
      const double w = static_cast<const JF2_t&>(*this).getValue(args...);
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF1_t&>(gradient) *= 1.0/w;
      static_cast<JF2_t&>(gradient)  = static_cast<const JF2_t&>(*this).getGradient(args...);
      static_cast<JF2_t&>(gradient) *= -v/(w*w);
 
      return gradient;
    }
  };
 
 
  /**
   * Auxiliary base class for fixed power of function.
   */
  template<>
  struct JFn<void>
  {};
 
  /**
   * Fixed power of function.
   */
  template<class JF1_t>
  struct JFn :
    public JFn<>,
    public JMathlib < JFn<JF1_t> >,
    public JF1_t
  { 
    using JMathlib< JFn<JF1_t> >::operator();
    using JMathlib< JFn<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JFn() :
      JF1_t(),
      N(1)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            first  function
     * \param  N             power
     */
    JFn(const JF1_t& f1, const int N) :
      JF1_t(f1),
      N(N)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(args...);
 
      return pow(u, N);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
      const double v = static_cast<const JF1_t&>(*this).getDerivative(x);
 
      return N * pow(u, N - 1) * v;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JFn getGradient(const Args& ...args) const
    {
      JFn gradient;
 
      const double u = static_cast<const JF1_t&>(*this).getValue(args...);
      const double w = N * pow(u, N - 1);
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(args...);
      static_cast<JF1_t&>(gradient) *= w;
 
      return gradient;
    }
 
  private:
    int N;
  };
 
 
  /**
   * x.
   */
  struct JX :
    public JMathlib < JX >,
    public JCalculus< JX >
  {
    /**
     * Default constructor.
     */
    JX()
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return x;
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return 1.0;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JX getGradient(const double x) const
    {
      JX gradient;
 
      return gradient;
    }
 
    static constexpr parameter_list<JX, 0> parameters{};
  };
 
 
  static const JX X;                                         //! Function object for x;
 
 
  /**
   * Fixed power of x.
   */
  template<int N>
  struct JXn :
    public JMathlib < JXn<N> >,
    public JCalculus< JXn<N> >
  {
    /**
     * Default constructor.
     */
    JXn()
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return pow(x, N);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return N * pow(x, N - 1);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JXn getGradient(const double x) const
    {
      JXn gradient;
 
      return gradient;
    }
 
    static constexpr parameter_list<JXn, 0> parameters{};
  };
 
 
  template<>
  struct JXn<0>;                                             //! Invalidate x^0.
  static const JXn<2> X2;                                    //! Function object for x^2;
  static const JXn<3> X3;                                    //! Function object for x^3;
  static const JXn<4> X4;                                    //! Function object for x^4;
  static const JXn<5> X5;                                    //! Function object for x^5;
 
 
  /**
   * Auxiliary base class for pointer to function.
   */
  template<>
  struct JPf<void>
  {};
 
  /**
   * Pointer to function.
   */
  template<class JF1_t>
  struct JPf :
    public JPf<>,
    public JMathlib < JPf<JF1_t> >,
    public JCalculus< JPf<JF1_t> >
  {
    using JMathlib< JPf<JF1_t> >::operator();
    using JMathlib< JPf<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JPf() :
      f1(NULL)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            pointer to function
     */
    JPf(const JF1_t* f1) :
      f1(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  args          abscissa value(s)
     * \return               function value
     */
    template<class ...Args>
    double getValue(const Args& ...args) const
    {
      return f1->getValue(args...);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return f1->getDerivative(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  args          abscissa value(s)
     * \return               gradient
     */
    template<class ...Args>
    JPf getGradient(const Args& ...args) const
    {
      JPf gradient;
 
      return gradient;
    }
 
    static constexpr parameter_list<JPf, 0> parameters{};
 
  private:
    const JF1_t* f1;
  };
 
 
  /**
   * Recursive template class for polynomial function.
   */
  template<int ID_t, size_t N>
  struct JPolynome :
    public JMathlib < JPolynome<ID_t, N> >,
    public JCalculus< JPolynome<ID_t, N> >,
    public JPolynome<ID_t, N - 1>
  {
    static const int ID = ID_t;
 
    static const size_t NUMBER_OF_DEGREES  =  N;
 
    using JCalculus< JPolynome<ID_t, N> >::negate;
    using JCalculus< JPolynome<ID_t, N> >::add;
    using JCalculus< JPolynome<ID_t, N> >::sub;
    using JCalculus< JPolynome<ID_t, N> >::mul;
    using JCalculus< JPolynome<ID_t, N> >::div;
    using JMathlib < JPolynome<ID_t, N> >::operator();
    using JMathlib < JPolynome<ID_t, N> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JPolynome() :
      a(0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JPolynome(const Args& ...args) :
      JPolynome()
    {
      set(args...);
    }
 
 
    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    JPolynome(const std::array<double, N+1>& args) :
      JPolynome()
    {
      for (size_t i = 0; i <= N; ++i) {
        (*this)[i] = std::data(args)[i];
      }
    }
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return  a * pow(x,N)  +  static_cast<const JPolynome<ID_t, N - 1>&>(*this).getValue(x); 
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return  N * a * pow(x,N-1)  +  static_cast<const JPolynome<ID_t, N - 1>&>(*this).getDerivative(x); 
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JPolynome getGradient(const double x) const
    {
      JPolynome gradient;
 
      gradient.a = pow(x,N);
      
      static_cast<JPolynome<ID_t, N - 1>&>(gradient) = static_cast<const JPolynome<ID_t, N - 1>&>(*this).getGradient(x);
 
      return gradient;
    }
 
    double a;                                                 //!< a[N]
 
    static constexpr parameter_list<JPolynome, N + 1> parameters = join<JPolynome>(JPolynome<ID_t, N-1>::parameters, &JPolynome::a);
 
  protected:
    /**
     * Recursive method for setting values.
     *
     * \param  x             first value
     * \param  args          remaining values
     */
    template<class ...Args>
    void set(const double x, const Args& ...args)
    {
      a = x;
 
      JPolynome<ID_t, N-1>::set(args...);
    }
  };
 
 
  /**
   * Termination class for polynomial function.
   * This function consititutes a constant term and is thereby multi-dimensional by nature.
   */
  template<int ID_t>
  struct JPolynome<ID_t, 0> :
    public JMathlib < JPolynome<ID_t, 0> >,
    public JCalculus< JPolynome<ID_t, 0> >
  {
    static const int ID = ID_t;
 
    static const size_t NUMBER_OF_DEGREES  =  0;
 
 
    /**
     * Default constructor.
     */
    JPolynome() :
      a(0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  a             value
     */
    JPolynome(const double a) :
      a(a)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    JPolynome(const std::array<double, 1>& args) :
      JPolynome()
    {
      a = std::data(args)[0];
    }
 
 
    /**
     * Function value.
     *
     * \return               function value
     */
    double getValue(...) const
    {
      return a;
    }
 
 
    /**
     * Derivative value.
     *
     * \return               derivative value
     */
    double getDerivative(...) const
    {
      return 0.0;
    }
 
 
    /**
     * Get gradient.
     *
     * \return               gradient
     */
    JPolynome getGradient(...) const
    {
      JPolynome gradient;
 
      gradient.a = 1.0;
 
      return gradient;
    }
 
    double a;                                                 //!< a[0]
 
    static constexpr parameter_list<JPolynome, 1> parameters{&JPolynome::a};
 
  protected:
    /**
     * Recursive method for setting values.
     *
     * \param  x             value
     */
    void set(const double x)
    {
      a = x;
    }
  };
 
 
  template<int ID_t>
  using JP0 = JPolynome<ID_t, 0>;                             //!< shorthand for 0th degree polynome
 
  template<int ID_t>
  using P0  = JP0<ID_t>;                                      //!< shorthand for constant term
 
  template<int ID_t>
  using JP1 = JPolynome<ID_t, 1>;                             //!< shorthand for 1st degree polynome
 
  template<int ID_t>
  using JP2 = JPolynome<ID_t, 2>;                             //!< shorthand for 2nd degree polynome
 
  template<int ID_t>
  using JP3 = JPolynome<ID_t, 3>;                             //!< shorthand for 3rd degree polynome
 
 
  /**
   * Gauss function.
   */
  template<int ID_t, bool normalised = false>
  struct JGauss :
    public JMathlib < JGauss<ID_t> >,
    public JCalculus< JGauss<ID_t> >
  {
    static const int ID = ID_t;
 
 
    /**
     * Default constructor.
     */
    JGauss() :
      center(0.0),
      sigma (0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  center          center
     * \param  sigma         sigma
     */
    JGauss(const double center,
           const double sigma) :
      center(center),
      sigma (sigma)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u  = (x - center) / sigma;
 
      return get(u);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double w  = 1.0 / sigma;
      const double u  = (x - center) / sigma;
 
      return get(u) * -u * w;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JGauss getGradient(const double x) const
    {
      JGauss gradient;
 
      const double w  = 1.0 / sigma;
      const double u  = (x - center) * w;
      const double f0 =  get(u);
 
      gradient.center =  f0 * u * w;                          // d(f)/d(center)
      gradient.sigma  =  f0 * u * u * w;                      // d(f)/d(sigma)
 
      return gradient;
    }
 
    double center;                                            //!< center
    double sigma;                                             //!< sigma
 
    static constexpr parameter_list<JGauss, 2> parameters{ &JGauss::center, &JGauss::sigma };
 
  private:
    /**
     * Get ordinate value.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      return exp(-0.5*u*u);
    }
  };
 
 
  /**
   * Gauss function.
   */
  template<int ID_t>
  struct JGauss<ID_t, true> :
    public JMathlib < JGauss<ID_t, true> >,
    public JCalculus< JGauss<ID_t, true> >
  {
    static const int ID = ID_t;
 
 
    /**
     * Default constructor.
     */
    JGauss() :
      center(0.0),
      sigma (0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  center          center
     * \param  sigma         sigma
     */
    JGauss(const double center,
           const double sigma) :
      center(center),
      sigma (sigma)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u  = (x - center) / sigma;
 
      return get(u);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double w  = 1.0 / sigma;
      const double u  = (x - center) / sigma;
 
      return get(u) * -u * w;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JGauss getGradient(const double x) const
    {
      JGauss gradient;
 
      const double w  = 1.0 / sigma;
      const double u  = (x - center) * w;
      const double f0 =  get(u);
 
      gradient.center =  f0 * (u)                   * w;      // d(f)/d(center)
      gradient.sigma  =  f0 * (u + 1.0) * (u - 1.0) * w;      // d(f)/d(sigma)
 
      return gradient;
    }
 
    double center;                                            //!< center
    double sigma;                                             //!< sigma
 
    static constexpr parameter_list<JGauss, 2> parameters{ &JGauss::center, &JGauss::sigma };
 
  private:
    /**
     * Get ordinate value.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      static const double W = 1.0 / sqrt(2.0 * acos(-1.0));
 
      return exp(-0.5*u*u) * W / sigma;
    }
  };
 
 
  /**
   * Power of function.
   */
  template<int ID_t, class JF1_t = JX>
  struct JPow :
    public JMathlib < JPow<ID_t, JF1_t> >,
    public JF1_t
  {
    static const int ID = ID_t;
 
    using JMathlib< JPow<ID_t, JF1_t> >::operator();
    using JMathlib< JPow<ID_t, JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JPow() :
      JF1_t(),
      alpha(0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  alpha         value
     */
    JPow(const double alpha) :
      JF1_t(),
      alpha(alpha)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     * \param  alpha         value
     */
    JPow(const JF1_t& f1, const double alpha) :
      JF1_t(f1),
      alpha(alpha)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      return pow(u, alpha);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
      const double v = static_cast<const JF1_t&>(*this).getDerivative(x);
 
      return alpha * pow(u, alpha - 1) * v;
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JPow getGradient(const double x) const
    {
      JPow gradient;
 
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
      const double w = pow(u, alpha);
 
      gradient.alpha = w * log(u);
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(x);
      static_cast<JF1_t&>(gradient) *= alpha * w / u;
 
      return gradient;
    }
 
    double alpha;                                             //!< f(x)^alpha
 
    static constexpr parameter_list<JPow, 1> parameters{ &JPow::alpha };
  };
 
 
  /**
   * Power of function.
   *
   * \param  f1            function
   * \param  alpha         alpha
   * \return               result
   */
  template<class JF1_t>
  JPow<JF1_t::ID, JF1_t> Pow(const JF1_t& f1, const double alpha)
  {
    return JPow<JF1_t::ID, JF1_t>(f1, alpha);
  }
 
 
  /**
   * Square root of function.
   */
  template<class JF1_t = JX>
  struct JSqrt :
    public JMathlib < JSqrt<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JSqrt<JF1_t> >::operator();
    using JMathlib< JSqrt<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JSqrt() :
      JF1_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JSqrt(const JF1_t& f1) :
      JF1_t(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      return sqrt(u);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
      const double v = static_cast<const JF1_t&>(*this).getDerivative(x);
 
      return 0.5 * v / sqrt(u);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JSqrt getGradient(const double x) const
    {
      JSqrt gradient;
 
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(x);
      static_cast<JF1_t&>(gradient) *= 0.5 / sqrt(u);
 
      return gradient;
    }
  };
 
 
  /**
   * Square root of function.
   *
   * \param  f1            function
   * \return               result
   */
  template<class JF1_t>
  JSqrt<JF1_t> Sqrt(const JF1_t& f1)
  {
    return JSqrt<JF1_t>(f1);
  }
 
 
  /**
   * Sine of function.
   */
  template<class JF1_t = JX>
  struct JSin :
    public JMathlib < JSin<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JSin<JF1_t> >::operator();
    using JMathlib< JSin<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JSin() :
      JF1_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JSin(const JF1_t& f1) :
      JF1_t(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      return sin(u);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
      const double v = static_cast<const JF1_t&>(*this).getDerivative(x);
 
      return v * cos(u);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JSin getGradient(const double x) const
    {
      JSin gradient;
 
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(x);
      static_cast<JF1_t&>(gradient) *= cos(u);
 
      return gradient;
    }
  };
 
 
  /**
   * Sine of function.
   *
   * \param  f1            function
   * \return               result
   */
  template<class JF1_t>
  JSin<JF1_t> Sin(const JF1_t& f1)
  {
    return JSin<JF1_t>(f1);
  }
 
 
  /**
   * Cosine of function.
   */
  template<class JF1_t = JX>
  struct JCos :
    public JMathlib < JCos<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JCos<JF1_t> >::operator();
    using JMathlib< JCos<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JCos() :
      JF1_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JCos(const JF1_t& f1) :
      JF1_t(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      return cos(u);
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
      const double v = static_cast<const JF1_t&>(*this).getDerivative(x);
 
      return v * -sin(u);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JCos getGradient(const double x) const
    {
      JCos gradient;
 
      const double u = static_cast<const JF1_t&>(*this).getValue(x);
 
      static_cast<JF1_t&>(gradient)  = static_cast<const JF1_t&>(*this).getGradient(x);
      static_cast<JF1_t&>(gradient) *= -sin(u);
 
      return gradient;
    }
  };
 
 
  /**
   * Cosine of function.
   *
   * \param  f1            function
   * \return               result
   */
  template<class JF1_t>
  JCos<JF1_t> Cos(const JF1_t& f1)
  {
    return JCos<JF1_t>(f1);
  }
 
 
  /**
   * Exponent of function.
   */
  template<class JF1_t = JX>
  struct JExp :
    public JMathlib< JExp<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JExp<JF1_t> >::operator();
    using JMathlib< JExp<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JExp() :
      JF1_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JExp(const JF1_t& f1) :
      JF1_t(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return exp(static_cast<const JF1_t&>(*this).getValue(x));
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return static_cast<const JF1_t&>(*this).getDerivative(x) * getValue(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JExp getGradient(const double x) const
    {
      JExp gradient;
 
      gradient  = static_cast<const JF1_t&>(*this).getGradient(x);
      gradient *= getValue(x);
 
      return gradient;
    }
  };
 
 
  /**
   * Exponent of zero degree polynomial function.
   */
  template<int ID_t>
  struct JExp< JPolynome<ID_t, 0> > :
    public JMathlib< JExp< JPolynome<ID_t, 0> > >,
    public JPolynome<ID_t, 0>
  {
    using JMathlib< JExp< JPolynome<ID_t, 0> > >::operator();
    using JMathlib< JExp< JPolynome<ID_t, 0> > >::operator[];
 
 
    /**
     * Default constructor.
     */
    JExp() :
      JPolynome<ID_t, 0>()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JExp(const JPolynome<ID_t, 0>& f1) :
      JPolynome<ID_t, 0>(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \return               function value
     */
    double getValue(...) const
    {
      return exp(static_cast<const JPolynome<ID_t, 0>&>(*this).getValue());
    }
 
 
    /**
     * Derivative value.
     *
     * \return               derivative value
     */
    double getDerivative(...) const
    {
      return static_cast<const JPolynome<ID_t, 0>&>(*this).getDerivative() * getValue();
    }
 
 
    /**
     * Get gradient.
     *
     * \return               gradient
     */
    JExp getGradient(...) const
    {
      JExp gradient;
 
      gradient  = static_cast<const JPolynome<ID_t, 0>&>(*this).getGradient();
      gradient *= getValue();
 
      return gradient;
    }
  };
 
 
  /**
   * Exponent of function.
   *
   * \param  f1            function
   * \return               result
   */
  template<class JF1_t>
  JExp<JF1_t> Exp(const JF1_t& f1)
  {
    return JExp<JF1_t>(f1);
  }
 
 
  /**
   * Logarithm of function.
   */
  template<class JF1_t = JX>
  struct JLog :
    public JMathlib< JLog<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JLog<JF1_t> >::operator();
    using JMathlib< JLog<JF1_t> >::operator[];
 
 
    /**
     * Default constructor.
     */
    JLog() :
      JF1_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JLog(const JF1_t& f1) :
      JF1_t(f1)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return log(static_cast<const JF1_t&>(*this).getValue(x));
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return static_cast<const JF1_t&>(*this).getDerivative(x) / static_cast<const JF1_t&>(*this).getValue(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    JLog getGradient(const double x) const
    {
      JLog gradient;
 
      gradient  =  static_cast<const JF1_t&>(*this).getGradient(x);
      gradient /= (static_cast<const JF1_t&>(*this).getValue(x));
 
      return gradient;
    }
  };
 
 
  /**
   * Logarithm of function.
   *
   * \param  f1            function
   * \return               result
   */
  template<class JF1_t>
  JLog<JF1_t> Log(const JF1_t& f1)
  {
    return JLog<JF1_t>(f1);
  }
}
 
#endif