#ifndef __JMATH__JMATHLIB__
#define __JMATH__JMATHLIB__

#include <cstddef>
#include <vector>
#include <initializer_list>
#include <cmath>

#include "JMath/JConstants.hh"


/**
 * \file
 * Functional algebra.
 *
 * \author mdejong
 */

namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  struct _vF;                //!< void function


  /**
   * Auxiliary data structure to make C-array from list of values.
   */
  struct make_carray :
    public std::vector<double>
  {
    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    make_carray(const Args& ...args)
    {
      add(args...);
    }
    

    /**
     * Type conversion operator.
     *
     * \return               pointer to list of values
     */
    operator const double*() const
    {
      return this->data();
    }

  private:
    /**
     * Recursive method for adding values to array.
     *
     * \param  x             first value
     * \param  args          remaining values
     */
    template<class ...Args>
    void add(const double x, const Args& ...args)
    {
      this->push_back(x);

      add(args...);
    }


    /**
     * Termination method for adding values to array.
     *
     * \param  x             final value
     */
    void add(const double x)
    {
      this->push_back(x);
    }
  };


  /**
   * Auxiliary data structure 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>
  struct parameter_list :
    public std::vector<double JF1_t::*>
  {
    /**
     * Default constructor.
     */
    parameter_list()
    {}


    /**
     * Constructor.
     *
     * \param  parameter     parameter
     */
    parameter_list(double JF1_t::*parameter) :
      std::vector<double JF1_t::*>(1, parameter)
    {}


    /**
     * Constructor.
     *
     * \param  parameters    parameters
     */
    parameter_list(const std::initializer_list<double JF1_t::*>& parameters) :
      std::vector<double JF1_t::*>(parameters)
    {}


    /**
     * Combine constructor.
     *
     * \param  parameters    parameters
     * \param  parameter     parameter
     */
    template<class U>
    parameter_list(const parameter_list<U>& parameters, double JF1_t::* parameter) :
      std::vector<double JF1_t::*>(parameters.begin(), parameters.end())
    {
      this->push_back(parameter);
    }


    /**
     * Combine constructor.
     *
     * \param  first         first  parameters
     * \param  second        second parameters
     */
    template<class U, class V>
    parameter_list(const parameter_list<U>& first, const parameter_list<V>& second)
    {
      for (const auto& i : first)  { this->push_back(i); }
      for (const auto& i : second) { this->push_back(i); }
    }
  };


  /**
   * Get number of parameters.
   *
   * \return               number of parameters
   */
  template<class JF1_t>
  inline size_t getNumberOfParameters()
  {
    return JF1_t::parameters.size();
  }


  /**
   * Get value of parameter at given index.
   *
   * \param  f1            function
   * \param  i             index
   * \return               value
   */
  template<class JF1_t>
  inline double getParameter(const JF1_t& f1, const size_t i)
  {
    return f1.*JF1_t::parameters.at(i);
  }


  /**
   * 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 != JF1_t::parameters.size(); ++i) {
	f1->*JF1_t::parameters[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 negate 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>
  struct JFn;                //!< forward declaration for fixed power of 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  function      this function
     * \param  N             power
     * \return               result function
     */
    friend JFn<JF1_t> operator^(const JF1_t& function, int N)
    {
      return JFn<JF1_t>(function, N);
    }
  };


  /**
   * Negate of function.
   */
  template<class JF1_t>
  struct JNegate :
    public JMathlib< JNegate<JF1_t> >,
    public JF1_t
  {
    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>
    const JNegate& getGradient(const Args& ...args) const
    {
      static 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();


    /**
     * 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>
    const JAdd& getGradient(const Args& ...args) const
    {
      static 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();


    /**
     * 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>
    const JSub& getGradient(const Args& ...args) const
    {
      static 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();


    /**
     * 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>
    const JMul& getGradient(const Args& ...args) const
    {
      static 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();


    /**
     * 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>
    const JDiv& getGradient(const Args& ...args) const
    {
      static 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>
  struct JPair :
    public JF1_t,
    public JF2_t,
    public JCalculus< JPair<JF1_t, JF2_t> >
  {
    using JCalculus< JPair<JF1_t, JF2_t> >::negate;
    using JCalculus< JPair<JF1_t, JF2_t> >::add;
    using JCalculus< JPair<JF1_t, JF2_t> >::sub;
    using JCalculus< JPair<JF1_t, JF2_t> >::mul;
    using JCalculus< JPair<JF1_t, JF2_t> >::div;


    /**
     * 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 const parameter_list<JPair> parameters;            //!< parameters
  };

  /**
   * Set parameters.
   */
  template<class JF1_t, class JF2_t>
  const parameter_list< JPair<JF1_t, JF2_t> > JPair<JF1_t, JF2_t>::parameters(JF1_t::parameters, JF2_t::parameters);


  /**
   * Addition of two functions.
   */
  template<class JF1_t, class JF2_t>
  struct JAdd :
    public JMathlib< JAdd<JF1_t, JF2_t> >,
    public JPair<JF1_t, JF2_t>
  {
    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<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>
    const JAdd& getGradient(const Args& ...args) const
    {
      static 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<JF1_t, JF2_t>
  {
    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<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>
    const JSub& getGradient(const Args& ...args) const
    {
      static 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<JF1_t, JF2_t>
  {
    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<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>
    const JMul& getGradient(const Args& ...args) const
    {
      static 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<JF1_t, JF2_t>
  {
    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<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>
    const JDiv& getGradient(const Args& ...args) const
    {
      static 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;
    }
  };


  /**
   * Fixed power of x.
   */
  template<class JF1_t>
  struct JFn :
    public JMathlib < JFn<JF1_t> >,
    public JF1_t
  { 

    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>
    const JFn& getGradient(const Args& ...args) const
    {
      static 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;
    }

    int N;
  };


  /**
   * 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();


    /**
     * Default constructor.
     */
    JPolynome() :
      a(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JPolynome(const Args& ...args) :
      JPolynome()
    {
      set(make_carray(args...));
    }


    /**
     * 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
     */
    const JPolynome& getGradient(const double x) const
    {
      static JPolynome gradient;                              // d(f)/d(a)

      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 const parameter_list<JPolynome> parameters;        //!< parameters

  protected:
    /**
     * Set parameter values.
     *
     * \param  pA            pointer to list of values
     */
    void set(const double* pA)
    {
      a = pA[N];

      JPolynome<ID_t, N-1>::set(pA);
    }
  };

  /**
   * Set parameters.
   */
  template<int ID_t, size_t N>
  const parameter_list< JPolynome<ID_t, N> > JPolynome<ID_t, N>::parameters(JPolynome<ID_t, N - 1>::parameters, &JPolynome<ID_t, N>::a);


  /**
   * Termination class for polynomial function.
   */
  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)
    {}


    /**
     * 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
     */
    const JPolynome& getGradient(...) const
    {
      static JPolynome gradient;                              // d(f)/d(a)

      gradient.a = 1.0;

      return gradient;
    }

    double a;                                                 //!< a[0]
    static const parameter_list<JPolynome> parameters;        //!< parameters

  protected:
    /**
     * Set parameter values.
     *
     * \param  array         pointer to list of values
     */
    void set(const double* array)
    {
      a = array[0];
    }
  };

  /**
   * Set parameters.
   */
  template<int ID_t>
  const parameter_list< JPolynome<ID_t, 0> > JPolynome<ID_t, 0>::parameters(&JPolynome<ID_t, 0>::a);


  template<int ID_t>
  using JP0 = JPolynome<ID_t, 0>;                             //!< short hand for 0th degree polynome

  template<int ID_t>
  using JP1 = JPolynome<ID_t, 1>;                             //!< short hand for 1st degree polynome

  template<int ID_t>
  using JP2 = JPolynome<ID_t, 2>;                             //!< short hand for 2nd degree polynome

  template<int ID_t>
  using JP3 = JPolynome<ID_t, 3>;                             //!< short hand 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
     */
    const JGauss& getGradient(const double x) const
    {
      static 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 const parameter_list<JGauss> parameters;           //!< parameters

  private:
    /**
     * Get ordinate value.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      return exp(-0.5*u*u);
    }
  };

  /**
   * Set parameters.
   */
  template<int ID_t, bool normalised>
  const parameter_list< JGauss<ID_t, normalised> > JGauss<ID_t, normalised>::parameters = { &JGauss<ID_t, normalised>::center, &JGauss<ID_t, normalised>::sigma };


  /**
   * 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
     */
    const JGauss& getGradient(const double x) const
    {
      static 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 const parameter_list<JGauss> parameters;           //!< parameters

  private:
    /**
     * Get ordinate value.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      return exp(-0.5*u*u) / (sqrt(2.0*PI) * sigma);
    }
  };

  /**
   * Set parameters.
   */
  template<int ID_t>
  const parameter_list< JGauss<ID_t, true> > JGauss<ID_t, true>::parameters = { &JGauss<ID_t, true>::center, &JGauss<ID_t, true>::sigma };


  /**
   * Power of function.
   */
  template<int ID_t, class JF1_t = _vF>
  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();


    /**
     * Default constructor.
     */
    JPow() :
      JF1_t(),
      alpha(0.0)
    {}


    /**
     * 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
     */
    const JPow& getGradient(const double x) const
    {
      static JPow gradient;                                   // d(f)/d(a)

      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 const parameter_list<JPow> parameters;             //!< parameters
  };

  /**
   * Set parameters.
   */
  template<int ID_t, class JF1_t>
  const parameter_list< JPow<ID_t, JF1_t> > JPow<ID_t, JF1_t>::parameters(JF1_t::parameters, &JPow<ID_t, JF1_t>::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);
  }


  /**
   * Power of x.
   */
  template<int ID_t>
  struct JPow<ID_t, _vF> :
    public JMathlib < JPow<ID_t> >,
    public JCalculus< JPow<ID_t> >
  {
    static const int ID = ID_t;


    /**
     * Default constructor.
     */
    JPow() :
      alpha(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  alpha         value
     */
    JPow(const double alpha) :
      alpha(alpha)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return pow(x, alpha);
    }


    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return alpha * pow(x, alpha - 1);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    const JPow& getGradient(const double x) const
    {
      static JPow gradient;                                   // d(f)/d(a)

      gradient.alpha = this->getValue(x) * log(x);

      return gradient;
    }

    double alpha;                                             //!< x^alpha
    static const parameter_list<JPow> parameters;             //!< parameters
  };

  /**
   * Set parameters.
   */
  template<int ID_t>
  const parameter_list< JPow<ID_t> > JPow<ID_t>::parameters(&JPow<ID_t>::alpha);


  /**
   * 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
     */
    const JXn& getGradient(const double x) const
    {
      static JXn gradient;

      return gradient;
    }

    static const parameter_list<JXn> parameters;              //!< parameters
  };

  /**
   * Set parameters.
   */
  template<int N>
  const parameter_list< JXn<N> > JXn<N>::parameters;


  /**
   * Fixed power of x.
   */
  template<>
  struct JXn<0>;


  /**
   * Square root of function.
   */
  template<class JF1_t = _vF>
  struct JSqrt :
    public JMathlib < JSqrt<JF1_t> >,
    public JF1_t
  {
    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
     */
    const JSqrt& getGradient(const double x) const
    {
      static JSqrt gradient;                                  // d(f)/d(a)

      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);
  }


  /**
   * Square root of x.
   */
  template<>
  struct JSqrt<_vF> :
    public JMathlib < JSqrt<_vF> >,
    public JCalculus< JSqrt<_vF> >
  {
    using JMathlib< JSqrt<_vF> >::operator();


    /**
     * Default constructor.
     */
    JSqrt()
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      return sqrt(x);
    }


    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      return 0.5 / sqrt(x);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    const JSqrt& getGradient(const double x) const
    {
      static JSqrt gradient;

      return gradient;
    }

    static const parameter_list<JSqrt> parameters;            //!< parameters
  };

  /**
   * Set parameters.
   */
  const parameter_list< JSqrt<> > JSqrt<>::parameters;


  /**
   * Square root of x.
   *
   * \return               result
   */
  JSqrt<> Sqrt()
  {
    return JSqrt<>();
  }


  /**
   * Sine of function.
   */
  template<class JF1_t>
  struct JSin :
    public JMathlib < JSin<JF1_t> >,
    public JF1_t
  {
    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
     */
    const JSin& getGradient(const double x) const
    {
      static JSin gradient;                                   // d(f)/d(a)

      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>
  struct JCos :
    public JMathlib < JCos<JF1_t> >,
    public JF1_t
  {
    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
     */
    const JCos& getGradient(const double x) const
    {
      static JCos gradient;                                   // d(f)/d(a)

      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 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>
  struct JExp :
    public JMathlib< JExp<JF1_t> >,
    public JF1_t
  {
    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
     */
    const JExp& getGradient(const double x) const
    {
      static 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();

    
    /**
     * 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
     */
    const JExp& getGradient(...) const
    {
      static 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>
  struct JLog :
    public JMathlib< JLog<JF1_t> >,
    public JF1_t
  {
    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
     */
    const JLog& getGradient(const double x) const
    {
      static 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