adtb_types.h

Go to the documentation of this file.

/*----------------------------------------------------------------------------
 ADOL-C -- Automatic Differentiation by Overloading in C++
 File:     adtb_types.h
 Revision: $Id$
 Contents: adtb_types.h contains the declarations of tape_location, adouble and
pdouble.
 
 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_AD_TYPES_H
#define ADOLC_AD_TYPES_H
 
#include <adolc/adolcerror.h>
#include <adolc/adolcexport.h>
#include <adolc/tape_interface.h>
#include <adolc/valuetape/valuetape.h>
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <span>
#include <type_traits>
 
// Forward declarations for concept
class adouble;
class pdouble;
 
template <typename T>
concept adouble_or_pdouble =
    std::is_same_v<T, adouble> || std::is_same_v<T, pdouble>;
 
template <adouble_or_pdouble T> class tape_location {
  size_t loc_{0}; 
 
  int valid_{1};
 
  size_t next_loc() {
    if constexpr (std::is_same_v<T, adouble>)
      return currentTape().next_loc();
    else
      return currentTape().p_next_loc();
  }
 
  void free_loc() {
    assert(currentTapePtr() != nullptr &&
           "Tape was deleted before all adoubles are destroyed!");
    if constexpr (std::is_same_v<T, adouble>)
      currentTape().free_loc(loc());
 
    // currently freeing the location of pdouble leads to errors
    // need to revisit the issue
    // currentTape().p_free_loc(loc());
  }
 
public:
  enum ownership_mode {
    INVALID = 0,
    OWNING = 1,
    BORROWED = 2,
  };
 
  ~tape_location() {
    if (owns_location()) {
      free_loc();
    }
  }
 
  tape_location() : loc_(next_loc()) {}
 
  explicit tape_location(size_t loc, ownership_mode mode = OWNING)
      : loc_(loc), valid_(mode) {}
 
  tape_location(const tape_location &) = delete;
 
  tape_location &operator=(const tape_location &) = delete;
 
  tape_location(tape_location &&other) noexcept
      : loc_(other.loc_), valid_(other.valid_) {
    other.valid_ = INVALID;
  };
 
  tape_location &operator=(tape_location &&other) noexcept {
    if (this == &other)
      return *this;
 
    // free location of *this
    if (owns_location())
      free_loc();
 
    // Transfer ownership
    loc_ = other.loc_;
    valid_ = other.valid_;
 
    other.valid_ = INVALID;
    return *this;
  };
 
  size_t loc() const { return loc_; }
 
  bool owns_location() const { return valid_ == OWNING; }
  bool is_valid() const { return valid_ != INVALID; }
 
  size_t release() {
    assert(is_valid());
    valid_ = INVALID;
    return loc_;
  }
};
 
class ADOLC_API adouble {
  tape_location<adouble> tape_loc_;
 
  explicit adouble(tape_location<adouble> &&loc) : tape_loc_(std::move(loc)) {}
 
public:
  ~adouble() = default;
 
  adouble();
 
  adouble(double coval);
 
  adouble(const adouble &a);
 
  adouble(adouble &&other) noexcept : tape_loc_(std::move(other.tape_loc_)) {};
 
  static adouble borrow_location(size_t loc) {
    return adouble(
        tape_location<adouble>(loc, tape_location<adouble>::BORROWED));
  }
 
  size_t release_loc() { return tape_loc_.release(); }
  bool is_valid() const { return tape_loc_.is_valid(); }
 
  // Assignment Operators
 
  adouble &operator=(const double coval);
 
  adouble &operator=(const adouble &a);
 
  adouble &operator=(adouble &&other) noexcept {
    if (this == &other) {
      return *this;
    }
    tape_loc_ = std::move(other.tape_loc_);
    return *this;
  }
  adouble &operator=(const pdouble &p);
 
  // Accessors
  double value() const {
    assert(is_valid());
    return currentTape().get_ad_value(loc());
  }
 
  size_t loc() const {
    assert(is_valid());
    return tape_loc_.loc();
  }
 
  // Mutators
  void value(const double coval) {
    assert(is_valid());
    currentTape().set_ad_value(loc(), coval);
  }
 
  // Type Conversions
 
  explicit operator double() const { return value(); }
 
  // Arithmetic assignment Operators
 
  adouble &operator+=(const double coval);
  adouble &operator+=(const adouble &a);
  adouble &operator+=(adouble &&a);
  inline adouble &operator+=(const pdouble &p);
 
  adouble &operator-=(const double coval);
  adouble &operator-=(const adouble &a);
  adouble &operator-=(adouble &&a);
  inline adouble &operator-=(const pdouble &p);
 
  adouble &operator*=(const double coval);
  adouble &operator*=(const adouble &a);
  inline adouble &operator*=(const pdouble &p);
 
  inline adouble &operator/=(const double coval);
  inline adouble &operator/=(const adouble &a);
  inline adouble &operator/=(adouble &&a);
  inline adouble &operator/=(const pdouble &p);
 
  // Postfix and Prefix Increment and Decrement Operators
 
  adouble operator++(int);
  adouble operator--(int);
 
  adouble &operator++();
  adouble &operator--();
 
  // Independent and Dependent Markers
 
  adouble &operator<<=(const double input);
 
  adouble &operator>>=(double &out);
 
  void declareIndependent();
 
  void declareDependent();
};
 
class ADOLC_API pdouble {
  tape_location<pdouble> tape_loc_;
 
public:
  ~pdouble() = default;
 
  pdouble() = delete;
 
  pdouble(const pdouble &) = delete;
  pdouble &operator=(const pdouble &) = delete;
 
  pdouble(pdouble &&other) noexcept : tape_loc_{std::move(other.tape_loc_)} {};
 
  pdouble &operator=(pdouble &&other) noexcept {
    if (this == &other)
      return *this;
 
    tape_loc_ = std::move(other.tape_loc_);
    return *this;
  }
 
  explicit pdouble(const double pval);
 
  explicit operator adouble() const;
 
  size_t loc() const { return tape_loc_.loc(); }
 
  double value() const { return currentTape().get_pd_value(loc()); }
 
  void value(const double pval) { currentTape().set_pd_value(loc(), pval); }
};
 
ADOLC_API inline void operator<<=(std::span<adouble> indeps,
                                  std::span<const double> inputs) {
  assert(indeps.size() == inputs.size() && "Size mismatch!");
  for (size_t i = 0; i < indeps.size(); ++i) {
    indeps[i] <<= inputs[i];
  }
}
 
ADOLC_API inline void operator>>=(std::span<adouble> deps,
                                  std::span<double> outputs) {
  assert(deps.size() == outputs.size() && "Size mismatch!");
  for (size_t i = 0; i < deps.size(); ++i) {
    deps[i] >>= outputs[i];
  }
}
 
ADOLC_API std::ostream &operator<<(std::ostream &, const adouble &);
ADOLC_API std::istream &operator>>(std::istream &, const adouble &);
 
/*--------------------------------------------------------------------------*/
/* sign operators */
ADOLC_API adouble operator+(const adouble &a);
ADOLC_API adouble operator+(adouble &&a);
 
ADOLC_API adouble operator-(const adouble &a);
ADOLC_API adouble operator-(adouble &&a);
 
/*--------------------------------------------------------------------------*/
/* binary operators */
 
/*--------------------------------------------------------------------------*/
ADOLC_API adouble operator+(const adouble &a, const adouble &b);
ADOLC_API adouble operator+(adouble &&a, const adouble &b);
ADOLC_API inline adouble operator+(const adouble &a, adouble &&b) {
  return std::move(b) + a;
}
ADOLC_API inline adouble operator+(adouble &&a, adouble &&b) {
  return std::move(a) + b;
}
ADOLC_API adouble operator+(const double coval, const adouble &a);
ADOLC_API adouble operator+(const double coval, adouble &&a);
inline adouble operator+(const adouble &a, const double coval) {
  return coval + a;
}
inline adouble operator+(adouble &&a, const double coval) {
  return coval + std::move(a);
}
/*--------------------------------------------------------------------------*/
 
ADOLC_API adouble operator-(const adouble &a, const adouble &b);
ADOLC_API adouble operator-(adouble &&a, const adouble &b);
ADOLC_API inline adouble operator-(const adouble &a, adouble &&b) {
  return -(std::move(b)) + a;
}
ADOLC_API inline adouble operator-(adouble &&a, adouble &&b) {
  return std::move(a) - b;
}
ADOLC_API adouble operator-(const double coval, const adouble &a);
ADOLC_API adouble operator-(const double coval, adouble &&a);
ADOLC_API inline adouble operator-(const adouble &a, const double coval) {
  return (-coval) + a;
}
ADOLC_API inline adouble operator-(adouble &&a, const double coval) {
  return (-coval) + std::move(a);
}
 
/*--------------------------------------------------------------------------*/
 
ADOLC_API adouble operator*(const adouble &a, const adouble &b);
ADOLC_API adouble operator*(adouble &&a, const adouble &b);
ADOLC_API inline adouble operator*(const adouble &a, adouble &&b) {
  return std::move(b) * a;
}
ADOLC_API inline adouble operator*(adouble &&a, adouble &&b) {
  return std::move(a) * b;
}
ADOLC_API adouble operator*(const double coval, const adouble &a);
ADOLC_API adouble operator*(const double coval, adouble &&a);
ADOLC_API inline adouble operator*(const adouble &a, const double coval) {
  return coval * a;
}
ADOLC_API inline adouble operator*(adouble &&a, const double coval) {
  return coval * std::move(a);
}
 
/*--------------------------------------------------------------------------*/
 
ADOLC_API adouble operator/(const adouble &a, const adouble &b);
ADOLC_API adouble operator/(adouble &&a, const adouble &b);
ADOLC_API adouble operator/(const adouble &a, adouble &&b);
ADOLC_API inline adouble operator/(adouble &&a, adouble &&b) {
  return std::move(a) / b;
}
ADOLC_API adouble operator/(const double coval, const adouble &a);
ADOLC_API adouble operator/(const double coval, adouble &&a);
ADOLC_API inline adouble operator/(const adouble &a, const double coval) {
  return a * (1.0 / coval);
}
ADOLC_API inline adouble operator/(adouble &&a, const double coval) {
  return std::move(a) * (1.0 / coval);
}
 
/*--------------------------------------------------------------------------*/
/* Inlined Arithmetic assignments  */
ADOLC_API inline adouble &adouble::operator/=(const double coval) {
  *this *= 1.0 / coval;
  return *this;
}
ADOLC_API inline adouble &adouble::operator/=(const adouble &a) {
  *this *= 1.0 / a;
  return *this;
}
ADOLC_API inline adouble &adouble::operator/=(adouble &&a) {
  *this *= 1.0 / std::move(a);
  return *this;
}
/*--------------------------------------------------------------------------*/
/* Comparison  */
 
ADOLC_API bool operator!=(const adouble &a, const double coval);
ADOLC_API bool operator!=(const double coval, const adouble &a);
 
ADOLC_API bool operator==(const double coval, const adouble &a);
ADOLC_API bool operator==(const adouble &a, const double coval);
 
ADOLC_API bool operator<=(const double coval, const adouble &a);
ADOLC_API bool operator<=(const adouble &a, const double coval);
 
ADOLC_API bool operator>=(const double coval, const adouble &a);
ADOLC_API bool operator>=(const adouble &a, const double coval);
 
ADOLC_API bool operator<(const double coval, const adouble &a);
ADOLC_API bool operator<(const adouble &a, const double coval);
 
ADOLC_API bool operator>(const double coval, const adouble &a);
ADOLC_API bool operator>(const adouble &a, const double coval);
 
#ifdef ADOLC_ADVANCED_BRANCHING
 
ADOLC_API adouble operator!=(const adouble &a, const adouble &b);
ADOLC_API adouble operator!=(adouble &&a, const adouble &b);
ADOLC_API adouble operator!=(const adouble &a, adouble &&b);
 
ADOLC_API adouble operator==(const adouble &a, const adouble &b);
ADOLC_API adouble operator==(adouble &&a, const adouble &b);
ADOLC_API adouble operator==(const adouble &a, adouble &&b);
 
ADOLC_API adouble operator<=(const adouble &a, const adouble &b);
ADOLC_API adouble operator<=(adouble &&a, const adouble &b);
ADOLC_API adouble operator<=(const adouble &a, adouble &&b);
 
ADOLC_API adouble operator>=(const adouble &a, const adouble &b);
ADOLC_API adouble operator>=(adouble &&a, const adouble &b);
ADOLC_API adouble operator>=(const adouble &a, adouble &&b);
 
ADOLC_API adouble operator<(const adouble &a, const adouble &b);
ADOLC_API adouble operator<(adouble &&a, const adouble &b);
ADOLC_API adouble operator<(const adouble &a, adouble &&b);
 
ADOLC_API adouble operator>(const adouble &a, const adouble &b);
ADOLC_API adouble operator>(adouble &&a, const adouble &b);
ADOLC_API adouble operator>(const adouble &a, adouble &&b);
 
#else // ADOLC_ADVANCED_BRANCHING
 
ADOLC_API inline bool operator!=(const adouble &a, const adouble &b) {
  return (a - b != 0.0);
}
ADOLC_API inline bool operator==(const adouble &a, const adouble &b) {
  return (a - b == 0.0);
}
ADOLC_API inline bool operator<=(const adouble &a, const adouble &b) {
  return (a - b <= 0.0);
}
ADOLC_API inline bool operator>=(const adouble &a, const adouble &b) {
  return (a - b >= 0.0);
}
ADOLC_API inline bool operator>(const adouble &a, const adouble &b) {
  return (a - b > 0.0);
}
ADOLC_API inline bool operator<(const adouble &a, const adouble &b) {
  return (a - b < 0.0);
}
 
#endif // ADOLC_ADVANCED_BRANCHING
 
ADOLC_API inline bool isnormal(const adouble &a) {
  return std::isnormal(a.value());
}
 
ADOLC_API inline bool isnan(const adouble &a) { return std::isnan(a.value()); }
 
ADOLC_API inline bool isinf(const adouble &a) { return std::isinf(a.value()); }
 
ADOLC_API inline bool isfinite(const adouble &a) {
  return std::isfinite(a.value());
}
 
/*--------------------------------------------------------------------------*/
/* unary operators */
ADOLC_API adouble exp(const adouble &a);
ADOLC_API adouble exp(adouble &&a);
 
ADOLC_API adouble log(const adouble &a);
ADOLC_API adouble log(adouble &&a);
 
ADOLC_API inline adouble log10(const adouble &a) {
  return log(a) / ADOLC_MATH_NSP::log(10.0);
}
ADOLC_API inline adouble log10(adouble &&a) {
  return log(std::move(a)) / ADOLC_MATH_NSP::log(10.0);
}
 
ADOLC_API adouble sqrt(const adouble &a);
ADOLC_API adouble sqrt(adouble &&a);
 
ADOLC_API adouble cbrt(const adouble &a);
ADOLC_API adouble cbrt(adouble &&a);
 
ADOLC_API adouble sin(const adouble &a);
ADOLC_API adouble sin(adouble &&a);
 
ADOLC_API adouble cos(const adouble &a);
ADOLC_API adouble cos(adouble &&a);
 
ADOLC_API inline adouble tan(const adouble &a) { return sin(a) / cos(a); }
ADOLC_API inline adouble tan(adouble &&a) { return sin(a) / cos(a); }
 
ADOLC_API adouble asin(const adouble &a);
ADOLC_API adouble asin(adouble &&a);
 
ADOLC_API adouble acos(const adouble &a);
ADOLC_API adouble acos(adouble &&a);
 
ADOLC_API adouble atan(const adouble &a);
ADOLC_API adouble atan(adouble &&a);
 
/*--------------------------------------------------------------------------*/
ADOLC_API inline adouble sinh(const adouble &a) {
  if (a.value() < 0.0) {
    adouble temp = exp(a);
    return 0.5 * (temp - 1.0 / temp);
  } else {
    adouble temp = exp(-a);
    return 0.5 * (1.0 / temp - temp);
  }
}
ADOLC_API inline adouble sinh(adouble &&a) {
  if (a.value() < 0.0) {
    adouble temp = exp(std::move(a));
    return 0.5 * (temp - 1.0 / temp);
  } else {
    adouble temp = exp(-std::move(a));
    return 0.5 * (1.0 / temp - temp);
  }
}
/*--------------------------------------------------------------------------*/
 
ADOLC_API inline adouble cosh(const adouble &a) {
  adouble temp = (a.value() < 0.0) ? exp(a) : exp(-a);
  return 0.5 * (temp + 1.0 / temp);
}
ADOLC_API inline adouble cosh(adouble &&a) {
  adouble temp = (a.value() < 0.0) ? exp(std::move(a)) : exp(-std::move(a));
  return 0.5 * (temp + 1.0 / temp);
}
/*--------------------------------------------------------------------------*/
 
ADOLC_API inline adouble tanh(const adouble &a) {
  if (a.value() < 0.0) {
    adouble temp = exp(2.0 * a);
    return (temp - 1.0) / (temp + 1.0);
  } else {
    adouble temp = exp((-2.0) * a);
    return (1.0 - temp) / (temp + 1.0);
  }
}
ADOLC_API inline adouble tanh(adouble &&a) {
  if (a.value() < 0.0) {
    adouble temp = exp(2.0 * std::move(a));
    return (temp - 1.0) / (temp + 1.0);
  } else {
    adouble temp = exp((-2.0) * std::move(a));
    return (1.0 - temp) / (temp + 1.0);
  }
}
/*--------------------------------------------------------------------------*/
 
ADOLC_API adouble asinh(const adouble &a);
ADOLC_API adouble asinh(adouble &&a);
 
ADOLC_API adouble acosh(const adouble &a);
ADOLC_API adouble acosh(adouble &&a);
 
ADOLC_API adouble atanh(const adouble &a);
ADOLC_API adouble atanh(adouble &&a);
 
ADOLC_API adouble erf(const adouble &a);
ADOLC_API adouble erf(adouble &&a);
 
ADOLC_API adouble erfc(const adouble &a);
ADOLC_API adouble erfc(adouble &&a);
 
ADOLC_API adouble ceil(const adouble &a);
ADOLC_API adouble ceil(adouble &&a);
 
ADOLC_API adouble floor(const adouble &a);
ADOLC_API adouble floor(adouble &&a);
 
ADOLC_API adouble fabs(const adouble &a);
ADOLC_API adouble fabs(adouble &&a);
 
ADOLC_API inline adouble abs(const adouble &a) { return fabs(a); }
ADOLC_API inline adouble abs(adouble &&a) { return fabs(std::move(a)); }
 
/*--------------------------------------------------------------------------*/
/* Min and Max */
 
ADOLC_API adouble fmin(const adouble &a, const adouble &b);
ADOLC_API adouble fmin(adouble &&a, const adouble &b);
ADOLC_API adouble fmin(const adouble &a, adouble &&b);
ADOLC_API inline adouble fmin(adouble &&a, adouble &&b) {
  return fmin(std::move(a), b);
}
ADOLC_API inline adouble fmin(const double coval, const adouble &a) {
  return fmin(adouble(coval), a);
}
ADOLC_API inline adouble fmin(const double coval, adouble &&a) {
  return fmin(adouble(coval), std::move(a));
}
ADOLC_API inline adouble fmin(const adouble &a, const double coval) {
  return fmin(a, adouble(coval));
}
ADOLC_API inline adouble fmin(adouble &&a, const double coval) {
  return fmin(std::move(a), adouble(coval));
}
 
ADOLC_API inline adouble min(const adouble &a, const adouble &b) {
  return fmin(a, b);
}
ADOLC_API inline adouble min(adouble &&a, const adouble &b) {
  return fmin(std::move(a), b);
}
ADOLC_API inline adouble min(const adouble &a, adouble &&b) {
  return fmin(a, std::move(b));
}
ADOLC_API inline adouble min(adouble &&a, adouble &&b) {
  return fmin(std::move(a), b);
}
ADOLC_API inline adouble min(const double coval, const adouble &a) {
  return fmin(coval, a);
}
ADOLC_API inline adouble min(const double coval, adouble &&a) {
  return fmin(coval, std::move(a));
}
ADOLC_API inline adouble min(const adouble &a, const double coval) {
  return fmin(a, coval);
}
ADOLC_API inline adouble min(adouble &&a, const double coval) {
  return fmin(std::move(a), coval);
}
 
/*--------------------------------------------------------------------------*/
 
ADOLC_API inline adouble fmax(const adouble &a, const adouble &b) {
  return (-fmin(-a, -b));
}
ADOLC_API inline adouble fmax(adouble &&a, const adouble &b) {
  return (-fmin(-std::move(a), -b));
}
ADOLC_API inline adouble fmax(const adouble &a, adouble &&b) {
  return (-fmin(-a, -std::move(b)));
}
ADOLC_API inline adouble fmax(adouble &&a, adouble &&b) {
  return fmax(std::move(a), b);
}
ADOLC_API inline adouble fmax(const double coval, const adouble &a) {
  return (-fmin(-coval, -a));
}
ADOLC_API inline adouble fmax(const double coval, adouble &&a) {
  return (-fmin(-coval, -std::move(a)));
}
ADOLC_API inline adouble fmax(const adouble &a, const double coval) {
  return (-fmin(-a, -coval));
}
ADOLC_API inline adouble fmax(adouble &&a, const double coval) {
  return (-fmin(-std::move(a), -coval));
}
 
ADOLC_API inline adouble max(const adouble &a, const adouble &b) {
  return fmax(a, b);
}
 
ADOLC_API inline adouble max(adouble &&a, const adouble &b) {
  return fmax(std::move(a), b);
}
ADOLC_API inline adouble max(const adouble &a, adouble &&b) {
  return fmax(a, std::move(b));
}
 
ADOLC_API inline adouble max(adouble &&a, adouble &&b) {
  return fmax(std::move(a), b);
}
ADOLC_API inline adouble max(const double coval, const adouble &a) {
  return fmax(coval, a);
}
ADOLC_API inline adouble max(const double coval, adouble &&a) {
  return fmax(coval, std::move(a));
}
ADOLC_API inline adouble max(const adouble &a, const double coval) {
  return fmax(a, coval);
}
ADOLC_API inline adouble max(adouble &&a, const double coval) {
  return fmax(std::move(a), coval);
}
 
/*--------------------------------------------------------------------------*/
/* ldexp and frexp */
 
ADOLC_API inline adouble ldexp(const adouble &a, const int exp) {
  return a * ldexp(1.0, exp);
}
ADOLC_API inline adouble ldexp(adouble &&a, const int exp) {
  return a * ldexp(1.0, exp);
}
 
/*--------------------------------------------------------------------------*/
 
ADOLC_API inline adouble frexp(const adouble &a, int *exp) {
  double coval = std::frexp(a.value(), exp);
  return adouble(coval);
}
ADOLC_API inline adouble frexp(adouble &&a, int *exp) {
  const double coval = std::frexp(a.value(), exp);
  a.value(coval);
  return a;
}
 
/*---------------------------------------------------------------*/
/* Conditional assignments */
 
ADOLC_API void condassign(adouble &res, const adouble &cond,
                          const adouble &arg1, const adouble &arg2);
ADOLC_API inline void condassign(double &res, const double &cond,
                                 const double &arg1, const double &arg2) {
  res = cond > 0 ? arg1 : arg2;
}
/*---------------------------------------------------------------*/
ADOLC_API void condassign(adouble &res, const adouble &cond,
                          const adouble &arg);
ADOLC_API inline void condassign(double &res, const double &cond,
                                 const double &arg) {
  res = cond > 0 ? arg : res;
}
/*---------------------------------------------------------------*/
ADOLC_API void condeqassign(adouble &res, const adouble &cond,
                            const adouble &arg1, const adouble &arg2);
ADOLC_API inline void condeqassign(double &res, const double &cond,
                                   const double &arg1, const double &arg2) {
  res = cond >= 0 ? arg1 : arg2;
}
/*---------------------------------------------------------------*/
ADOLC_API inline void condeqassign(double &res, const double &cond,
                                   const double &arg) {
  res = cond >= 0 ? arg : res;
}
ADOLC_API void condeqassign(adouble &res, const adouble &cond,
                            const adouble &arg);
 
/*---------------------------------------------------------------*/
/* power function */
ADOLC_API adouble pow(const adouble &a, const double exp);
ADOLC_API adouble pow(adouble &&a, const double exp);
ADOLC_API inline adouble pow(double coval, const adouble &a) {
  adouble ret;
  if (coval <= 0)
    ADOLCError::fail(ADOLCError::ErrorType::NONPOSITIVE_BASIS,
                     CURRENT_LOCATION);
 
  condassign(ret, adouble{coval}, exp(a * ADOLC_MATH_NSP::log(coval)),
             adouble{ADOLC_MATH_NSP::pow(coval, a.value())});
 
  return ret;
}
ADOLC_API inline adouble pow(const adouble &a, const adouble &b) {
  assert((a.value() >= 0) && "\nADOL-C message: negative basis deactivated\n ");
  assert(a.value() != 0 && "\nADOL-C message: zero basis deactivated\n ");
 
  adouble a1, a2, ret;
 
  condassign(a1, -b, adouble{ADOLC_MATH_NSP::pow(a.value(), b.value())},
             pow(a, b.value()));
  condassign(a2, fabs(a), pow(a, b.value()), a1);
  condassign(ret, a, exp(b * log(a)), a2);
  return ret;
}
 
/*---------------------------------------------------------------*/
/* atan(x, y) */
ADOLC_API inline adouble atan2(const adouble &a, const adouble &b) {
  adouble a1, a2, ret, sy;
  const double pihalf = ADOLC_MATH_NSP::asin(1.0);
  condassign(sy, a, adouble{1.0}, adouble{-1.0});
  condassign(a1, b, atan(a / b), atan(a / b) + sy * 2 * pihalf);
  condassign(a2, fabs(a), sy * pihalf - atan(b / a), adouble{0.0});
  condassign(ret, fabs(b) - fabs(a), a1, a2);
  return ret;
}
 
/*---------------------------------------------------------------*/
/* Advector functions  */
// copies src to dest; both should be of size "size"
void ADOLC_API adolc_vec_copy(adouble *const dest, const adouble *const src,
                              size_t size);
 
// computes canonical scalarproduct, src and dest must have size "size"
void ADOLC_API adolc_vec_dot(adouble *const dest, const adouble *const src,
                             size_t size);
 
// res = a*vec_a + vec_b, vec_a and vec_b must have same size "size"
void ADOLC_API adolc_vec_axpy(adouble *const res, const adouble &a,
                              const adouble *const vec_a,
                              const adouble *const vec_b, size_t size);
 
/*--------------------------------------------------------------------------*/
/* Arithmetic assignments */
ADOLC_API inline adouble &adouble::operator+=(const pdouble &p) {
  return *this += adouble(p);
}
ADOLC_API inline adouble &adouble::operator-=(const pdouble &p) {
  return *this -= adouble(p);
}
ADOLC_API inline adouble &adouble::operator*=(const pdouble &p) {
  return *this *= adouble(p);
}
ADOLC_API inline adouble &adouble::operator/=(const pdouble &p) {
  return *this /= adouble(p);
}
 
/*--------------------------------------------------------------------------*/
/* sign operators */
ADOLC_API inline adouble operator-(const pdouble &p) { return -adouble(p); }
ADOLC_API inline adouble operator+(const pdouble &p) { return adouble(p); };
 
/*--------------------------------------------------------------------------*/
/* Binary Operations */
ADOLC_API inline adouble operator+(const pdouble &p, const adouble &a) {
  return a + adouble(p);
};
ADOLC_API inline adouble operator+(const adouble &a, const pdouble &p) {
  return a + adouble(p);
};
ADOLC_API inline adouble operator+(const pdouble &p, const double coval) {
  return coval + adouble(p);
}
ADOLC_API inline adouble operator+(const double coval, const pdouble &p) {
  return coval + adouble(p);
}
/*--------------------------------------------------------------------------*/
ADOLC_API inline adouble operator-(const adouble &a, const pdouble &p) {
  return a - adouble(p);
}
ADOLC_API inline adouble operator-(const pdouble &p, const adouble &a) {
  return (-a) + adouble(p);
}
ADOLC_API inline adouble operator-(const pdouble &p, const double coval) {
  return adouble(p) - coval;
}
ADOLC_API inline adouble operator-(const double coval, const pdouble &p) {
  return coval - adouble(p);
}
/*--------------------------------------------------------------------------*/
ADOLC_API inline adouble operator*(const adouble &a, const pdouble &p) {
  return a * adouble(p);
}
ADOLC_API inline adouble operator*(const pdouble &p, const adouble &a) {
  return a * adouble(p);
}
ADOLC_API inline adouble operator*(const pdouble &p, const double coval) {
  return coval * adouble(p);
}
ADOLC_API inline adouble operator*(const double coval, const pdouble &p) {
  return coval * adouble(p);
}
/*--------------------------------------------------------------------------*/
ADOLC_API inline adouble operator/(const pdouble &p, const adouble &a) {
  return adouble(p) / a;
}
ADOLC_API inline adouble operator/(const adouble &a, const pdouble &p) {
  return a / adouble(p);
}
ADOLC_API inline adouble operator/(const pdouble &p, const double coval) {
  return adouble(p) / coval;
}
ADOLC_API inline adouble operator/(const double coval, const pdouble &p) {
  return coval / adouble(p);
}
 
/*--------------------------------------------------------------------------*/
/* Comparions */
#ifdef ADOLC_ADVANCED_BRANCHING
ADOLC_API inline adouble operator!=(const adouble &a, const pdouble &p) {
  return a != adouble(p);
}
ADOLC_API inline adouble operator!=(const pdouble &p, const adouble &a) {
  return adouble(p) != a;
}
 
ADOLC_API inline adouble operator==(const adouble &a, const pdouble &p) {
  return a == adouble(p);
}
ADOLC_API inline adouble operator==(const pdouble &p, const adouble &a) {
  return adouble(p) == a;
}
 
ADOLC_API inline adouble operator<=(const adouble &a, const pdouble &p) {
  return a <= adouble(p);
}
ADOLC_API inline adouble operator<=(const pdouble &p, const adouble &a) {
  return adouble(p) <= a;
}
 
ADOLC_API inline adouble operator>=(const adouble &a, const pdouble &p) {
  return a >= adouble(p);
}
ADOLC_API inline adouble operator>=(const pdouble &p, const adouble &a) {
  return adouble(p) >= a;
}
 
ADOLC_API inline adouble operator>(const adouble &a, const pdouble &p) {
  return a > adouble(p);
}
ADOLC_API inline adouble operator>(const pdouble &p, const adouble &a) {
  return adouble(p) > a;
}
ADOLC_API inline adouble operator<(const adouble &a, const pdouble &p) {
  return a < adouble(p);
}
ADOLC_API inline adouble operator<(const pdouble &p, const adouble &a) {
  return adouble(p) < a;
}
#else // ADOLC_ADVANCED_BRANCHING
 
ADOLC_API inline bool operator!=(const adouble &a, const pdouble &p) {
  return a != adouble(p);
}
ADOLC_API inline bool operator!=(const pdouble &p, const adouble &a) {
  return adouble(p) != a;
}
 
ADOLC_API inline bool operator==(const adouble &a, const pdouble &p) {
  return a == adouble(p);
}
ADOLC_API inline bool operator==(const pdouble &p, const adouble &a) {
  return adouble(p) == a;
}
 
ADOLC_API inline bool operator<=(const adouble &a, const pdouble &p) {
  return a <= adouble(p);
}
ADOLC_API inline bool operator<=(const pdouble &p, const adouble &a) {
  return adouble(p) <= a;
}
 
ADOLC_API inline bool operator>=(const adouble &a, const pdouble &p) {
  return a >= adouble(p);
}
ADOLC_API inline bool operator>=(const pdouble &p, const adouble &a) {
  return adouble(p) >= a;
}
 
ADOLC_API inline bool operator>(const adouble &a, const pdouble &p) {
  return a > adouble(p);
}
ADOLC_API inline bool operator>(const pdouble &p, const adouble &a) {
  return adouble(p) > a;
}
 
ADOLC_API inline bool operator<(const adouble &a, const pdouble &p) {
  return a < adouble(p);
}
ADOLC_API inline bool operator<(const pdouble &p, const adouble &a) {
  return adouble(p) < a;
}
 
#endif // ADOLC_ADVANCED_BRANCHING
 
ADOLC_API inline bool isnormal(const pdouble &p) {
  return std::isnormal(p.value());
}
 
ADOLC_API inline bool isnan(const pdouble &p) { return std::isnan(p.value()); }
 
ADOLC_API inline bool isinf(const pdouble &p) { return std::isinf(p.value()); }
 
ADOLC_API inline bool isfinite(const pdouble &p) {
  return std::isfinite(p.value());
}
 
/*--------------------------------------------------------------------------*/
/* unary operators */
ADOLC_API inline adouble exp(const pdouble &p) { return exp(adouble(p)); }
ADOLC_API inline adouble log(const pdouble &p) { return log(adouble(p)); }
ADOLC_API inline adouble log10(const pdouble &p) { return log10(adouble(p)); }
ADOLC_API inline adouble sqrt(const pdouble &p) { return sqrt(adouble(p)); }
ADOLC_API inline adouble cbrt(const pdouble &p) { return cbrt(adouble(p)); }
ADOLC_API inline adouble sin(const pdouble &p) { return sin(adouble(p)); }
ADOLC_API inline adouble cos(const pdouble &p) { return cos(adouble(p)); }
ADOLC_API inline adouble tan(const pdouble &p) { return tan(adouble(p)); }
ADOLC_API inline adouble asin(const pdouble &p) { return asin(adouble(p)); }
ADOLC_API inline adouble acos(const pdouble &p) { return acos(adouble(p)); }
ADOLC_API inline adouble atan(const pdouble &p) { return atan(adouble(p)); }
ADOLC_API inline adouble sinh(const pdouble &p) { return sinh(adouble(p)); }
ADOLC_API inline adouble cosh(const pdouble &p) { return cosh(adouble(p)); }
ADOLC_API inline adouble tanh(const pdouble &p) { return tanh(adouble(p)); }
ADOLC_API inline adouble asinh(const pdouble &p) { return asinh(adouble(p)); }
ADOLC_API inline adouble acosh(const pdouble &p) { return acosh(adouble(p)); }
ADOLC_API inline adouble atanh(const pdouble &p) { return atanh(adouble(p)); }
ADOLC_API inline adouble erf(const pdouble &p) { return erf(adouble(p)); }
ADOLC_API inline adouble erfc(const pdouble &p) { return erfc(adouble(p)); }
ADOLC_API inline adouble fabs(const pdouble &p) { return fabs(adouble(p)); }
ADOLC_API inline adouble abs(const pdouble &p) { return fabs(adouble(p)); }
ADOLC_API inline adouble ceil(const pdouble &p) { return ceil(adouble(p)); }
ADOLC_API inline adouble floor(const pdouble &p) { return floor(adouble(p)); }
 
/*--------------------------------------------------------------------------*/
/* Min and Max */
ADOLC_API inline adouble fmin(const adouble &a, const pdouble &p) {
  return fmin(a, adouble(p));
}
ADOLC_API inline adouble fmin(const pdouble &p, const adouble &a) {
  return fmin(adouble(p), a);
}
ADOLC_API inline adouble fmin(const pdouble &p, const pdouble &q) {
  return fmin(adouble(p), adouble(q));
}
ADOLC_API inline adouble fmin(const double coval, const pdouble &p) {
  return fmin(coval, adouble(p));
}
ADOLC_API inline adouble fmin(const pdouble &p, const double coval) {
  return fmin(adouble(p), coval);
}
 
ADOLC_API inline adouble min(const adouble &a, const pdouble &p) {
  return fmin(a, p);
}
ADOLC_API inline adouble min(const pdouble &p, const adouble &a) {
  return fmin(p, a);
}
ADOLC_API inline adouble min(const pdouble &p, const pdouble &q) {
  return fmin(p, q);
}
ADOLC_API inline adouble min(const double coval, const pdouble &p) {
  return fmin(coval, p);
}
ADOLC_API inline adouble min(const pdouble &p, const double coval) {
  return fmin(p, coval);
}
 
/*--------------------------------------------------------------------------*/
ADOLC_API inline adouble fmax(const adouble &a, const pdouble &p) {
  return fmax(a, adouble(p));
}
ADOLC_API inline adouble fmax(const pdouble &p, const adouble &a) {
  return fmax(adouble(p), a);
}
ADOLC_API inline adouble fmax(const pdouble &p, const pdouble &q) {
  return fmax(adouble(p), adouble(q));
}
ADOLC_API inline adouble fmax(const double coval, const pdouble &p) {
  return fmax(coval, adouble(p));
}
ADOLC_API inline adouble fmax(const pdouble &p, const double coval) {
  return fmax(adouble(p), coval);
}
 
ADOLC_API inline adouble max(const adouble &a, const pdouble &p) {
  return fmax(a, p);
}
ADOLC_API inline adouble max(const pdouble &p, const adouble &a) {
  return fmax(p, a);
}
ADOLC_API inline adouble max(const pdouble &p, const pdouble &q) {
  return fmax(p, q);
}
ADOLC_API inline adouble max(const double coval, const pdouble &p) {
  return fmax(coval, p);
}
ADOLC_API inline adouble max(const pdouble &p, const double coval) {
  return fmax(p, coval);
}
 
/*--------------------------------------------------------------------------*/
/* ldexp and frexp */
ADOLC_API inline adouble ldexp(const pdouble &p, const int exp) {
  return ldexp(adouble(p), exp);
}
ADOLC_API inline adouble frexp(const pdouble &p, int *exp) {
  return frexp(adouble(p), exp);
}
/*--------------------------------------------------------------------------*/
/* power function */
ADOLC_API inline adouble pow(const adouble &a, const pdouble &p) {
  return pow(a, adouble(p));
}
ADOLC_API inline adouble pow(const pdouble &p, const adouble &a) {
  return pow(adouble(p), a);
}
ADOLC_API inline adouble pow(const pdouble &p, const double coval) {
  return pow(adouble(p), coval);
}
 
/*--------------------------------------------------------------------------*/
/* atan(x, y) */
ADOLC_API inline adouble atan2(const pdouble &p, const adouble &a) {
  return atan2(adouble(p), a);
}
ADOLC_API inline adouble atan2(const adouble &a, const pdouble &p) {
  return atan2(a, adouble(p));
}
 
//-----------------------------------------------------------
/* User defined version of logarithm to test extend_quad macro */
ADOLC_API double myquad(double x);
 
template <> struct ADOLC_API std::numeric_limits<adouble> {
  static constexpr bool is_specialized = true;
 
  static constexpr double min() noexcept {
    return std::numeric_limits<double>::min();
  }
 
  static constexpr double max() noexcept {
    return std::numeric_limits<double>::max();
  }
 
  static constexpr double lowest() noexcept {
    return std::numeric_limits<double>::lowest();
  }
 
  static constexpr int digits = std::numeric_limits<double>::digits;
  static constexpr int digits10 = std::numeric_limits<double>::digits10;
  static constexpr int max_digits10 = std::numeric_limits<double>::max_digits10;
  static constexpr bool is_signed = std::numeric_limits<double>::is_signed;
  static constexpr bool is_integer = std::numeric_limits<double>::is_integer;
  static constexpr bool is_exact = std::numeric_limits<double>::is_exact;
  static constexpr int radix = std::numeric_limits<double>::radix;
 
  static constexpr double epsilon() noexcept {
    return std::numeric_limits<double>::epsilon();
  }
 
  static constexpr double round_error() noexcept {
    return std::numeric_limits<double>::round_error();
  }
 
  static constexpr int min_exponent = std::numeric_limits<double>::min_exponent;
  static constexpr int min_exponent10 =
      std::numeric_limits<double>::min_exponent10;
  static constexpr int max_exponent = std::numeric_limits<double>::max_exponent;
  static constexpr int max_exponent10 =
      std::numeric_limits<double>::max_exponent10;
 
  static constexpr bool has_infinity =
      std::numeric_limits<double>::has_infinity;
  static constexpr bool has_quiet_NaN =
      std::numeric_limits<double>::has_quiet_NaN;
  static constexpr bool has_signaling_NaN =
      std::numeric_limits<double>::has_signaling_NaN;
  static constexpr float_denorm_style has_denorm =
      std::numeric_limits<double>::has_denorm;
  static constexpr bool has_denorm_loss =
      std::numeric_limits<double>::has_denorm_loss;
 
  static constexpr double infinity() noexcept {
    return std::numeric_limits<double>::infinity();
  }
 
  static constexpr double quiet_NaN() noexcept {
    return std::numeric_limits<double>::quiet_NaN();
  }
 
  static constexpr double signaling_NaN() noexcept {
    return std::numeric_limits<double>::signaling_NaN();
  }
 
  static constexpr double denorm_min() noexcept {
    return std::numeric_limits<double>::denorm_min();
  }
 
  static constexpr bool is_iec559 = std::numeric_limits<double>::is_iec559;
  static constexpr bool is_bounded = std::numeric_limits<double>::is_bounded;
  static constexpr bool is_modulo = std::numeric_limits<double>::is_modulo;
 
  static constexpr bool traps = std::numeric_limits<double>::traps;
  static constexpr bool tinyness_before =
      std::numeric_limits<double>::tinyness_before;
  static constexpr float_round_style round_style =
      std::numeric_limits<double>::round_style;
};
 
#endif // ADOLC_AD_TYPES_H