adtl.h

Go to the documentation of this file.

/*----------------------------------------------------------------------------
 ADOL-C -- Automatic Differentiation by Overloading in C++
 File:     adouble.cpp
 Revision: $Id$
 Contents: adtl.h contains that declaratins of procedures used to
           define various traceless adouble operations.
 
 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel,
               Benjamin Letschert, Kshitij Kulshreshtha
 
 This file is part of ADOL-C. This software is provided as open source.
 Any use, reproduction, or distribution of the software constitutes
 recipient's acceptance of the terms of the accompanying license file.
 
----------------------------------------------------------------------------*/
#ifndef ADOLC_ADTL_H
#define ADOLC_ADTL_H
 
#include <adolc/adolcexport.h>
#include <adolc/internal/common.h>
#include <list>
#include <ostream>
#include <stdexcept>
 
#if USE_BOOST_POOL
#include <boost/pool/pool_alloc.hpp>
#endif
 
using std::istream;
using std::list;
using std::logic_error;
using std::ostream;
 
namespace adtl {
 
ADOLC_API double makeNaN();
ADOLC_API double makeInf();
 
#ifdef USE_ADTL_REFCOUNTING
class adouble;
 
class refcounter {
private:
  static size_t refcnt;
  ADOLC_API friend void setNumDir(const size_t p);
  ADOLC_API friend class adouble;
 
public:
  refcounter() { ++refcnt; }
  ~refcounter() { --refcnt; }
  ADOLC_API inline static size_t getNumLiveVar() { return refcnt; }
};
#endif
 
// class func_ad {
// public:
//     virtual int operator() (int n, adouble *x, int m, adouble *y) = 0;
// };
 
class adouble {
public:
  inline adouble();
  inline adouble(const double v);
  inline adouble(const double v, const double *adv);
  inline adouble(const adouble &a);
  inline ~adouble();
 
  // sign
  inline adouble operator-() const;
  inline adouble operator+() const;
 
  // addition
  inline adouble operator+(const double v) const;
  inline adouble operator+(const adouble &a) const;
  ADOLC_API inline friend adouble operator+(const double v, const adouble &a);
 
  // subtraction
  inline adouble operator-(const double v) const;
  inline adouble operator-(const adouble &a) const;
  ADOLC_API inline friend adouble operator-(const double v, const adouble &a);
 
  // multiplication
  inline adouble operator*(const double v) const;
  inline adouble operator*(const adouble &a) const;
  ADOLC_API inline friend adouble operator*(const double v, const adouble &a);
 
  // division
  inline adouble operator/(const double v) const;
  inline adouble operator/(const adouble &a) const;
  ADOLC_API inline friend adouble operator/(const double v, const adouble &a);
 
  // inc/dec
  inline adouble operator++();
  inline adouble operator++(int);
  inline adouble operator--();
  inline adouble operator--(int);
 
  // functions
  ADOLC_API inline friend adouble tan(const adouble &a);
  ADOLC_API inline friend adouble exp(const adouble &a);
  ADOLC_API inline friend adouble log(const adouble &a);
  ADOLC_API inline friend adouble sqrt(const adouble &a);
  ADOLC_API inline friend adouble cbrt(const adouble &a);
  ADOLC_API inline friend adouble sin(const adouble &a);
  ADOLC_API inline friend adouble cos(const adouble &a);
  ADOLC_API inline friend adouble asin(const adouble &a);
  ADOLC_API inline friend adouble acos(const adouble &a);
  ADOLC_API inline friend adouble atan(const adouble &a);
 
  ADOLC_API inline friend adouble atan2(const adouble &a, const adouble &b);
  ADOLC_API inline friend adouble pow(const adouble &a, double v);
  ADOLC_API inline friend adouble pow(const adouble &a, const adouble &b);
  ADOLC_API inline friend adouble pow(double v, const adouble &a);
  ADOLC_API inline friend adouble log10(const adouble &a);
 
  ADOLC_API inline friend adouble sinh(const adouble &a);
  ADOLC_API inline friend adouble cosh(const adouble &a);
  ADOLC_API inline friend adouble tanh(const adouble &a);
  ADOLC_API inline friend adouble asinh(const adouble &a);
  ADOLC_API inline friend adouble acosh(const adouble &a);
  ADOLC_API inline friend adouble atanh(const adouble &a);
  ADOLC_API inline friend adouble fabs(const adouble &a);
  ADOLC_API inline friend adouble ceil(const adouble &a);
  ADOLC_API inline friend adouble floor(const adouble &a);
  ADOLC_API inline friend adouble fmax(const adouble &a, const adouble &b);
  ADOLC_API inline friend adouble fmax(double v, const adouble &a);
  ADOLC_API inline friend adouble fmax(const adouble &a, double v);
  ADOLC_API inline friend adouble fmin(const adouble &a, const adouble &b);
  ADOLC_API inline friend adouble fmin(double v, const adouble &a);
  ADOLC_API inline friend adouble fmin(const adouble &a, double v);
  ADOLC_API inline friend adouble ldexp(const adouble &a, const adouble &b);
  ADOLC_API inline friend adouble ldexp(const adouble &a, const double v);
  ADOLC_API inline friend adouble ldexp(const double v, const adouble &a);
  ADOLC_API inline friend double frexp(const adouble &a, int *v);
  ADOLC_API inline friend adouble erf(const adouble &a);
  ADOLC_API inline friend adouble erfc(const adouble &a);
 
  ADOLC_API inline friend void condassign(adouble &res, const adouble &cond,
                                          const adouble &arg1,
                                          const adouble &arg2);
  ADOLC_API inline friend void condassign(adouble &res, const adouble &cond,
                                          const adouble &arg);
  ADOLC_API inline friend void condeqassign(adouble &res, const adouble &cond,
                                            const adouble &arg1,
                                            const adouble &arg2);
  ADOLC_API inline friend void condeqassign(adouble &res, const adouble &cond,
                                            const adouble &arg);
 
  /*******************  nontemporary results  ***************************/
  // assignment
  inline adouble &operator=(const double v);
  inline adouble &operator=(const adouble &a);
 
  // addition
  inline adouble &operator+=(const double v);
  inline adouble &operator+=(const adouble &a);
 
  // subtraction
  inline adouble &operator-=(const double v);
  inline adouble &operator-=(const adouble &a);
 
  // multiplication
  inline adouble &operator*=(const double v);
  inline adouble &operator*=(const adouble &a);
 
  // division
  inline adouble &operator/=(const double v);
  inline adouble &operator/=(const adouble &a);
 
  // not
  inline bool operator!() const;
 
  // comparison
  inline bool operator!=(const adouble &) const;
  inline bool operator!=(const double) const;
  ADOLC_API inline friend bool operator!=(const double, const adouble &);
 
  inline bool operator==(const adouble &) const;
  inline bool operator==(const double) const;
  ADOLC_API inline friend bool operator==(const double, const adouble &);
 
  inline bool operator<=(const adouble &) const;
  inline bool operator<=(const double) const;
  ADOLC_API inline friend bool operator<=(const double, const adouble &);
 
  inline bool operator>=(const adouble &) const;
  inline bool operator>=(const double) const;
  ADOLC_API inline friend bool operator>=(const double, const adouble &);
 
  inline bool operator>(const adouble &) const;
  inline bool operator>(const double) const;
  ADOLC_API inline friend bool operator>(const double, const adouble &);
 
  inline bool operator<(const adouble &) const;
  inline bool operator<(const double) const;
  ADOLC_API inline friend bool operator<(const double, const adouble &);
 
  /*******************  getter / setter  ********************************/
  inline double getValue() const;
  inline void setValue(const double v);
  inline double *getADValue() const;
  inline void setADValue(const double *v);
 
  inline double getADValue(size_t p) const;
  inline void setADValue(size_t p, double v);
  inline explicit operator double const &() const;
  inline explicit operator double &&();
  inline explicit operator double();
 
public:
  /*******************  i/o operations  *********************************/
  friend ostream &operator<<(ostream &, const adouble &);
  friend istream &operator>>(istream &, adouble &);
 
private:
#if USE_BOOST_POOL
  ADOLC_API static boost::pool<boost::default_user_allocator_new_delete>
      *advalpool;
#endif
  double *adval;
#ifdef USE_ADTL_REFCOUNTING
  refcounter __rcnt;
#endif
  ADOLC_API inline static size_t numDir = 1;
  ADOLC_API inline friend void setNumDir(const size_t p);
  ADOLC_API inline friend size_t getNumDir();
};
 
} // namespace adtl
 
#include <cmath>
#include <iostream>
#include <limits>
 
namespace adtl {
 
ADOLC_API inline void setNumDir(const size_t p) {
#ifdef USE_ADTL_REFCOUNTING
  if (refcounter::refcnt > 0) {
    fprintf(DIAG_OUT,
            "ADOL-C Warning: Traceless: Setting numDir will not change the "
            "number of\n directional derivative in existing adoubles and may "
            "lead to erroneous results\n or memory corruption\n Number of "
            "currently existing adoubles = %zu\n",
            refcounter::refcnt);
  }
#endif
 
  if (p < 1) {
    fprintf(DIAG_OUT, "ADOL-C Error: Traceless: p < 1 not possible\n");
    abort();
  }
  adouble::numDir = p;
#if USE_BOOST_POOL
  if (adouble::advalpool != NULL) {
    delete adouble::advalpool;
    adouble::advalpool = NULL;
  }
  adouble::advalpool =
      new boost::pool<boost::default_user_allocator_new_delete>(
          (adouble::numDir + 1) * sizeof(double));
#endif
}
 
ADOLC_API inline size_t getNumDir() { return adouble::numDir; }
 
ADOLC_API inline double makeNaN() {
  return ADOLC_MATH_NSP::numeric_limits<double>::quiet_NaN();
}
 
ADOLC_API inline double makeInf() {
  return ADOLC_MATH_NSP::numeric_limits<double>::infinity();
}
 
#define FOR_I_EQ_0_LTEQ_NUMDIR for (size_t _i = 0; _i <= adouble::numDir; ++_i)
#define FOR_I_EQ_1_LTEQ_NUMDIR for (size_t _i = 1; _i <= adouble::numDir; ++_i)
#define ADVAL_I adval[_i]
#define PRIMAL_VALUE adval[0]
 
/*******************************  ctors  ************************************/
inline adouble::adouble() : adval(NULL) {
#if USE_BOOST_POOL
  adval = reinterpret_cast<double *>(advalpool->malloc());
#else
  adval = new double[adouble::numDir + 1];
#endif
  PRIMAL_VALUE = 0.;
}
 
inline adouble::adouble(const double v) : adval(NULL) {
#if USE_BOOST_POOL
  adval = reinterpret_cast<double *>(advalpool->malloc());
#else
  adval = new double[adouble::numDir + 1];
#endif
  PRIMAL_VALUE = v;
  FOR_I_EQ_1_LTEQ_NUMDIR
  ADVAL_I = 0.0;
}
 
inline adouble::adouble(const double v, const double *adv) : adval(NULL) {
#if USE_BOOST_POOL
  adval = reinterpret_cast<double *>(advalpool->malloc());
#else
  adval = new double[adouble::numDir + 1];
#endif
  PRIMAL_VALUE = v;
  FOR_I_EQ_1_LTEQ_NUMDIR
  ADVAL_I = adv[_i - 1];
}
 
inline adouble::adouble(const adouble &a) : adval(NULL) {
#if USE_BOOST_POOL
  adval = reinterpret_cast<double *>(advalpool->malloc());
#else
  adval = new double[adouble::numDir + 1];
#endif
  FOR_I_EQ_0_LTEQ_NUMDIR
  ADVAL_I = a.ADVAL_I;
}
 
/*******************************  dtors  ************************************/
inline adouble::~adouble() {
  if (adval != NULL)
#if USE_BOOST_POOL
    advalpool->free(adval);
#else
    delete[] adval;
#endif
}
 
/*************************  temporary results  ******************************/
// sign
inline adouble adouble::operator-() const {
  adouble tmp;
  FOR_I_EQ_0_LTEQ_NUMDIR
  tmp.ADVAL_I = -ADVAL_I;
  return tmp;
}
 
inline adouble adouble::operator+() const { return *this; }
 
// addition
inline adouble adouble::operator+(const double v) const {
  adouble tmp(PRIMAL_VALUE + v, adval + 1);
  return tmp;
}
 
inline adouble adouble::operator+(const adouble &a) const {
  adouble tmp;
  FOR_I_EQ_0_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I + a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble operator+(const double v, const adouble &a) {
  adouble tmp(v + a.PRIMAL_VALUE, a.adval + 1);
  return tmp;
}
 
// subtraction
inline adouble adouble::operator-(const double v) const {
  adouble tmp(PRIMAL_VALUE - v, adval + 1);
  return tmp;
}
 
inline adouble adouble::operator-(const adouble &a) const {
  adouble tmp;
  FOR_I_EQ_0_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I - a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble operator-(const double v, const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = v - a.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = -a.ADVAL_I;
  return tmp;
}
 
// multiplication
inline adouble adouble::operator*(const double v) const {
  adouble tmp;
  FOR_I_EQ_0_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I * v;
  return tmp;
}
 
inline adouble adouble::operator*(const adouble &a) const {
  adouble tmp;
  tmp.PRIMAL_VALUE = PRIMAL_VALUE * a.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I * a.PRIMAL_VALUE + PRIMAL_VALUE * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble operator*(const double v, const adouble &a) {
  adouble tmp;
  FOR_I_EQ_0_LTEQ_NUMDIR
  tmp.ADVAL_I = v * a.ADVAL_I;
  return tmp;
}
 
// division
inline adouble adouble::operator/(const double v) const {
  adouble tmp;
  FOR_I_EQ_0_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I / v;
  return tmp;
}
 
inline adouble adouble::operator/(const adouble &a) const {
  adouble tmp;
  tmp.PRIMAL_VALUE = PRIMAL_VALUE / a.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = (ADVAL_I * a.PRIMAL_VALUE - PRIMAL_VALUE * a.ADVAL_I) /
                (a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  return tmp;
}
 
ADOLC_API inline adouble operator/(const double v, const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = v / a.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = (-v * a.ADVAL_I) / (a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  return tmp;
}
 
// inc/dec
inline adouble adouble::operator++() {
  ++PRIMAL_VALUE;
  return *this;
}
 
inline adouble adouble::operator++(int) {
  adouble tmp;
  tmp.PRIMAL_VALUE = PRIMAL_VALUE++;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I;
  return tmp;
}
 
inline adouble adouble::operator--() {
  --PRIMAL_VALUE;
  return *this;
}
 
inline adouble adouble::operator--(int) {
  adouble tmp;
  tmp.PRIMAL_VALUE = PRIMAL_VALUE--;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = ADVAL_I;
  return tmp;
}
 
// functions
ADOLC_API inline adouble tan(const adouble &a) {
  adouble tmp;
  double tmp2;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::tan(a.PRIMAL_VALUE);
  tmp2 = ADOLC_MATH_NSP::cos(a.PRIMAL_VALUE);
  tmp2 *= tmp2;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble exp(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::exp(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp.PRIMAL_VALUE * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble log(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::log(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  if (a.PRIMAL_VALUE > 0)
    tmp.ADVAL_I = a.ADVAL_I / a.PRIMAL_VALUE;
  else if (a.PRIMAL_VALUE == 0 && a.ADVAL_I != 0.0) {
    int sign = (a.ADVAL_I < 0) ? -1 : 1;
    tmp.ADVAL_I = sign * makeInf();
  } else
    tmp.ADVAL_I = makeNaN();
  return tmp;
}
 
ADOLC_API inline adouble sqrt(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::sqrt(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  if (a.PRIMAL_VALUE > 0)
    tmp.ADVAL_I = a.ADVAL_I / (tmp.PRIMAL_VALUE * 2);
  else if (a.PRIMAL_VALUE == 0.0 && a.ADVAL_I != 0.0) {
    int sign = (a.ADVAL_I < 0) ? -1 : 1;
    tmp.ADVAL_I = sign * makeInf();
  } else
    tmp.ADVAL_I = makeNaN();
  return tmp;
}
 
ADOLC_API inline adouble cbrt(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::cbrt(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  if (a.PRIMAL_VALUE != 0.0)
    tmp.ADVAL_I = a.ADVAL_I / (tmp.PRIMAL_VALUE * tmp.PRIMAL_VALUE * 3);
  else if (a.PRIMAL_VALUE == 0.0 && a.ADVAL_I != 0.0) {
    int sign = (a.ADVAL_I < 0) ? -1 : 1;
    tmp.ADVAL_I = sign * makeInf();
  } else
    tmp.ADVAL_I = makeNaN();
  return tmp;
}
 
ADOLC_API inline adouble sin(const adouble &a) {
  adouble tmp;
  double tmp2;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::sin(a.PRIMAL_VALUE);
  tmp2 = ADOLC_MATH_NSP::cos(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble cos(const adouble &a) {
  adouble tmp;
  double tmp2;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::cos(a.PRIMAL_VALUE);
  tmp2 = -ADOLC_MATH_NSP::sin(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble asin(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::asin(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::sqrt(1 - a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble acos(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::acos(a.PRIMAL_VALUE);
  double tmp2 = -ADOLC_MATH_NSP::sqrt(1 - a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble atan(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::atan(a.PRIMAL_VALUE);
  double tmp2 = 1 + a.PRIMAL_VALUE * a.PRIMAL_VALUE;
  tmp2 = 1 / tmp2;
  if (tmp2 != 0)
    FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I * tmp2;
  else FOR_I_EQ_1_LTEQ_NUMDIR tmp.ADVAL_I = 0.0;
  return tmp;
}
 
ADOLC_API inline adouble atan2(const adouble &a, const adouble &b) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::atan2(a.PRIMAL_VALUE, b.PRIMAL_VALUE);
  double tmp2 = a.PRIMAL_VALUE * a.PRIMAL_VALUE;
  double tmp3 = b.PRIMAL_VALUE * b.PRIMAL_VALUE;
  double tmp4 = tmp3 / (tmp2 + tmp3);
  if (tmp4 != 0)
    FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I =
      (a.ADVAL_I * b.PRIMAL_VALUE - a.PRIMAL_VALUE * b.ADVAL_I) / tmp3 * tmp4;
  else FOR_I_EQ_1_LTEQ_NUMDIR tmp.ADVAL_I = 0.0;
  return tmp;
}
 
ADOLC_API inline adouble pow(const adouble &a, double v) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::pow(a.PRIMAL_VALUE, v);
  double tmp2 = v * ADOLC_MATH_NSP::pow(a.PRIMAL_VALUE, v - 1);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble pow(const adouble &a, const adouble &b) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::pow(a.PRIMAL_VALUE, b.PRIMAL_VALUE);
  double tmp2 =
      b.PRIMAL_VALUE * ADOLC_MATH_NSP::pow(a.PRIMAL_VALUE, b.PRIMAL_VALUE - 1);
  double tmp3 = ADOLC_MATH_NSP::log(a.PRIMAL_VALUE) * tmp.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I + tmp3 * b.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble pow(double v, const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::pow(v, a.PRIMAL_VALUE);
  double tmp2 = tmp.PRIMAL_VALUE * ADOLC_MATH_NSP::log(v);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble log10(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::log10(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::log((double)10) * a.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble sinh(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::sinh(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::cosh(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I * tmp2;
  return tmp;
}
 
ADOLC_API inline adouble cosh(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::cosh(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::sinh(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I * tmp2;
  return tmp;
}
 
ADOLC_API inline adouble tanh(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::tanh(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::cosh(a.PRIMAL_VALUE);
  tmp2 *= tmp2;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble asinh(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP_ERF::asinh(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::sqrt(a.PRIMAL_VALUE * a.PRIMAL_VALUE + 1);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble acosh(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP_ERF::acosh(a.PRIMAL_VALUE);
  double tmp2 = ADOLC_MATH_NSP::sqrt(a.PRIMAL_VALUE * a.PRIMAL_VALUE - 1);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble atanh(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP_ERF::atanh(a.PRIMAL_VALUE);
  double tmp2 = 1 - a.PRIMAL_VALUE * a.PRIMAL_VALUE;
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I / tmp2;
  return tmp;
}
 
ADOLC_API inline adouble fabs(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::fabs(a.PRIMAL_VALUE);
  int as = 0;
  if (a.PRIMAL_VALUE > 0)
    as = 1;
  if (a.PRIMAL_VALUE < 0)
    as = -1;
  if (as != 0)
    FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = a.ADVAL_I * as;
  else FOR_I_EQ_1_LTEQ_NUMDIR {
    as = 0;
    if (a.ADVAL_I > 0)
      as = 1;
    if (a.ADVAL_I < 0)
      as = -1;
    tmp.ADVAL_I = a.ADVAL_I * as;
  }
  return tmp;
}
 
ADOLC_API inline adouble ceil(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::ceil(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = 0.0;
  return tmp;
}
 
ADOLC_API inline adouble floor(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP::floor(a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = 0.0;
  return tmp;
}
 
ADOLC_API inline adouble fmax(const adouble &a, const adouble &b) {
  adouble tmp;
  double tmp2 = a.PRIMAL_VALUE - b.PRIMAL_VALUE;
  if (tmp2 < 0) {
    tmp.PRIMAL_VALUE = b.PRIMAL_VALUE;
    FOR_I_EQ_1_LTEQ_NUMDIR
    tmp.ADVAL_I = b.ADVAL_I;
  } else {
    tmp.PRIMAL_VALUE = a.PRIMAL_VALUE;
    if (tmp2 > 0) {
      FOR_I_EQ_1_LTEQ_NUMDIR
      tmp.ADVAL_I = a.ADVAL_I;
    } else {
      FOR_I_EQ_1_LTEQ_NUMDIR {
        if (a.ADVAL_I < b.ADVAL_I)
          tmp.ADVAL_I = b.ADVAL_I;
        else
          tmp.ADVAL_I = a.ADVAL_I;
      }
    }
  }
  return tmp;
}
 
ADOLC_API inline adouble fmax(double v, const adouble &a) {
  adouble tmp;
  double tmp2 = v - a.PRIMAL_VALUE;
  if (tmp2 < 0) {
    tmp.PRIMAL_VALUE = a.PRIMAL_VALUE;
    FOR_I_EQ_1_LTEQ_NUMDIR
    tmp.ADVAL_I = a.ADVAL_I;
  } else {
    tmp.PRIMAL_VALUE = v;
    if (tmp2 > 0) {
      FOR_I_EQ_1_LTEQ_NUMDIR
      tmp.ADVAL_I = 0.0;
    } else {
      FOR_I_EQ_1_LTEQ_NUMDIR {
        if (a.ADVAL_I > 0)
          tmp.ADVAL_I = a.ADVAL_I;
        else
          tmp.ADVAL_I = 0.0;
      }
    }
  }
  return tmp;
}
 
ADOLC_API inline adouble fmax(const adouble &a, double v) {
  adouble tmp;
  double tmp2 = a.PRIMAL_VALUE - v;
  if (tmp2 < 0) {
    tmp.PRIMAL_VALUE = v;
    FOR_I_EQ_1_LTEQ_NUMDIR
    tmp.ADVAL_I = 0.0;
  } else {
    tmp.PRIMAL_VALUE = a.PRIMAL_VALUE;
    if (tmp2 > 0) {
      FOR_I_EQ_1_LTEQ_NUMDIR
      tmp.ADVAL_I = a.ADVAL_I;
    } else {
      FOR_I_EQ_1_LTEQ_NUMDIR {
        if (a.ADVAL_I > 0)
          tmp.ADVAL_I = a.ADVAL_I;
        else
          tmp.ADVAL_I = 0.0;
      }
    }
  }
  return tmp;
}
 
ADOLC_API inline adouble fmin(const adouble &a, const adouble &b) {
  adouble tmp;
  double tmp2 = a.PRIMAL_VALUE - b.PRIMAL_VALUE;
  if (tmp2 < 0) {
    tmp.PRIMAL_VALUE = a.PRIMAL_VALUE;
    FOR_I_EQ_1_LTEQ_NUMDIR
    tmp.ADVAL_I = a.ADVAL_I;
  } else {
    tmp.PRIMAL_VALUE = b.PRIMAL_VALUE;
    if (tmp2 > 0) {
      FOR_I_EQ_1_LTEQ_NUMDIR
      tmp.ADVAL_I = b.ADVAL_I;
    } else {
      FOR_I_EQ_1_LTEQ_NUMDIR {
        if (a.ADVAL_I < b.ADVAL_I)
          tmp.ADVAL_I = a.ADVAL_I;
        else
          tmp.ADVAL_I = b.ADVAL_I;
      }
    }
  }
  return tmp;
}
 
ADOLC_API inline adouble fmin(double v, const adouble &a) {
  adouble tmp;
  double tmp2 = v - a.PRIMAL_VALUE;
  if (tmp2 < 0) {
    tmp.PRIMAL_VALUE = v;
    FOR_I_EQ_1_LTEQ_NUMDIR
    tmp.ADVAL_I = 0.0;
  } else {
    tmp.PRIMAL_VALUE = a.PRIMAL_VALUE;
    if (tmp2 > 0) {
      FOR_I_EQ_1_LTEQ_NUMDIR
      tmp.ADVAL_I = a.ADVAL_I;
    } else {
      FOR_I_EQ_1_LTEQ_NUMDIR {
        if (a.ADVAL_I < 0)
          tmp.ADVAL_I = a.ADVAL_I;
        else
          tmp.ADVAL_I = 0.0;
      }
    }
  }
  return tmp;
}
 
ADOLC_API inline adouble fmin(const adouble &a, double v) {
  adouble tmp;
  double tmp2 = a.PRIMAL_VALUE - v;
  if (tmp2 < 0) {
    tmp.PRIMAL_VALUE = a.PRIMAL_VALUE;
    FOR_I_EQ_1_LTEQ_NUMDIR
    tmp.ADVAL_I = a.ADVAL_I;
  } else {
    tmp.PRIMAL_VALUE = v;
    if (tmp2 > 0) {
      FOR_I_EQ_1_LTEQ_NUMDIR
      tmp.ADVAL_I = 0.0;
    } else {
      FOR_I_EQ_1_LTEQ_NUMDIR {
        if (a.ADVAL_I < 0)
          tmp.ADVAL_I = a.ADVAL_I;
        else
          tmp.ADVAL_I = 0.0;
      }
    }
  }
  return tmp;
}
 
ADOLC_API inline adouble ldexp(const adouble &a, const adouble &b) {
  adouble tmp = a * pow(2., b);
  return tmp;
}
 
ADOLC_API inline adouble ldexp(const adouble &a, const double v) {
  return a * ADOLC_MATH_NSP::pow(2., v);
}
 
ADOLC_API inline adouble ldexp(const double v, const adouble &a) {
  adouble tmp = v * pow(2., a);
  return tmp;
}
 
ADOLC_API inline double frexp(const adouble &a, int *v) {
  return ADOLC_MATH_NSP::frexp(a.PRIMAL_VALUE, v);
}
 
ADOLC_API inline adouble erf(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP_ERF::erf(a.PRIMAL_VALUE);
  double tmp2 = 2.0 / ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
                ADOLC_MATH_NSP_ERF::exp(-a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline adouble erfc(const adouble &a) {
  adouble tmp;
  tmp.PRIMAL_VALUE = ADOLC_MATH_NSP_ERF::erfc(a.PRIMAL_VALUE);
  double tmp2 = -2.0 / ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
                ADOLC_MATH_NSP_ERF::exp(-a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  FOR_I_EQ_1_LTEQ_NUMDIR
  tmp.ADVAL_I = tmp2 * a.ADVAL_I;
  return tmp;
}
 
ADOLC_API inline void condassign(adouble &res, const adouble &cond,
                                 const adouble &arg1, const adouble &arg2) {
  if (cond.getValue() > 0)
    res = arg1;
  else
    res = arg2;
}
 
ADOLC_API inline void condassign(adouble &res, const adouble &cond,
                                 const adouble &arg) {
  if (cond.getValue() > 0)
    res = arg;
}
 
ADOLC_API inline void condeqassign(adouble &res, const adouble &cond,
                                   const adouble &arg1, const adouble &arg2) {
  if (cond.getValue() >= 0)
    res = arg1;
  else
    res = arg2;
}
 
ADOLC_API inline void condeqassign(adouble &res, const adouble &cond,
                                   const adouble &arg) {
  if (cond.getValue() >= 0)
    res = arg;
}
 
/*******************  nontemporary results  *********************************/
inline adouble &adouble::operator=(const double v) {
  PRIMAL_VALUE = v;
  FOR_I_EQ_1_LTEQ_NUMDIR
  ADVAL_I = 0.0;
  return *this;
}
 
inline adouble &adouble::operator=(const adouble &a) {
  FOR_I_EQ_0_LTEQ_NUMDIR
  ADVAL_I = a.ADVAL_I;
  return *this;
}
 
inline adouble &adouble::operator+=(const double v) {
  PRIMAL_VALUE += v;
  return *this;
}
 
inline adouble &adouble::operator+=(const adouble &a) {
  FOR_I_EQ_0_LTEQ_NUMDIR
  ADVAL_I += a.ADVAL_I;
  return *this;
}
 
inline adouble &adouble::operator-=(const double v) {
  PRIMAL_VALUE -= v;
  return *this;
}
 
inline adouble &adouble::operator-=(const adouble &a) {
  FOR_I_EQ_0_LTEQ_NUMDIR
  ADVAL_I -= a.ADVAL_I;
  return *this;
}
 
inline adouble &adouble::operator*=(const double v) {
  FOR_I_EQ_0_LTEQ_NUMDIR
  ADVAL_I *= v;
  return *this;
}
 
inline adouble &adouble::operator*=(const adouble &a) {
  FOR_I_EQ_1_LTEQ_NUMDIR
  ADVAL_I = ADVAL_I * a.PRIMAL_VALUE + PRIMAL_VALUE * a.ADVAL_I;
  PRIMAL_VALUE *= a.PRIMAL_VALUE;
  return *this;
}
 
inline adouble &adouble::operator/=(const double v) {
  FOR_I_EQ_0_LTEQ_NUMDIR
  ADVAL_I /= v;
  return *this;
}
 
inline adouble &adouble::operator/=(const adouble &a) {
  FOR_I_EQ_1_LTEQ_NUMDIR
  ADVAL_I = (ADVAL_I * a.PRIMAL_VALUE - PRIMAL_VALUE * a.ADVAL_I) /
            (a.PRIMAL_VALUE * a.PRIMAL_VALUE);
  PRIMAL_VALUE = PRIMAL_VALUE / a.PRIMAL_VALUE;
  return *this;
}
 
// not
inline bool adouble::operator!() const { return PRIMAL_VALUE == 0.0; }
 
// comparison
inline bool adouble::operator!=(const adouble &a) const {
  return PRIMAL_VALUE != a.PRIMAL_VALUE;
}
 
inline bool adouble::operator!=(const double v) const {
  return PRIMAL_VALUE != v;
}
 
ADOLC_API inline bool operator!=(const double v, const adouble &a) {
  return v != a.PRIMAL_VALUE;
}
 
inline bool adouble::operator==(const adouble &a) const {
  return PRIMAL_VALUE == a.PRIMAL_VALUE;
}
 
inline bool adouble::operator==(const double v) const {
  return PRIMAL_VALUE == v;
}
 
ADOLC_API inline bool operator==(const double v, const adouble &a) {
  return v == a.PRIMAL_VALUE;
}
 
inline bool adouble::operator<=(const adouble &a) const {
  return PRIMAL_VALUE <= a.PRIMAL_VALUE;
}
 
inline bool adouble::operator<=(const double v) const {
  return PRIMAL_VALUE <= v;
}
 
ADOLC_API inline bool operator<=(const double v, const adouble &a) {
  return v <= a.PRIMAL_VALUE;
}
 
inline bool adouble::operator>=(const adouble &a) const {
  return PRIMAL_VALUE >= a.PRIMAL_VALUE;
}
 
inline bool adouble::operator>=(const double v) const {
  return PRIMAL_VALUE >= v;
}
 
ADOLC_API inline bool operator>=(const double v, const adouble &a) {
  return v >= a.PRIMAL_VALUE;
}
 
inline bool adouble::operator>(const adouble &a) const {
  return PRIMAL_VALUE > a.PRIMAL_VALUE;
}
 
inline bool adouble::operator>(const double v) const {
  return PRIMAL_VALUE > v;
}
 
ADOLC_API inline bool operator>(const double v, const adouble &a) {
  return v > a.PRIMAL_VALUE;
}
 
inline bool adouble::operator<(const adouble &a) const {
  return PRIMAL_VALUE < a.PRIMAL_VALUE;
}
 
inline bool adouble::operator<(const double v) const {
  return PRIMAL_VALUE < v;
}
 
ADOLC_API inline bool operator<(const double v, const adouble &a) {
  return v < a.PRIMAL_VALUE;
}
 
/*******************  getter / setter  **************************************/
inline adouble::operator double const &() const { return PRIMAL_VALUE; }
 
inline adouble::operator double &&() { return (double &&)PRIMAL_VALUE; }
 
inline adouble::operator double() { return PRIMAL_VALUE; }
 
inline double adouble::getValue() const { return PRIMAL_VALUE; }
 
inline void adouble::setValue(const double v) { PRIMAL_VALUE = v; }
 
inline double *adouble::getADValue() const { return (adval + 1); }
 
inline void adouble::setADValue(const double *const v) {
  FOR_I_EQ_1_LTEQ_NUMDIR
  ADVAL_I = v[_i - 1];
}
 
inline double adouble::getADValue(size_t p) const {
  if (p >= adouble::numDir) {
    fprintf(DIAG_OUT, "Derivative array accessed out of bounds"
                      " while \"getADValue(...)\"!!!\n");
    throw logic_error("incorrect function call, errorcode=-1");
  }
  return adval[p + 1];
}
 
inline void adouble::setADValue(size_t p, double v) {
  if (p >= adouble::numDir) {
    fprintf(DIAG_OUT, "Derivative array accessed out of bounds"
                      " while \"setADValue(...)\"!!!\n");
    throw logic_error("incorrect function call, errorcode=-1");
  }
  adval[p + 1] = v;
}
 
} // namespace adtl
#endif