1866 lines
67 KiB
C++
1866 lines
67 KiB
C++
//===-- duals/dual - Dual number class --------------------------*- C++ -*-===//
|
||
//
|
||
// Part of the cppduals project.
|
||
// https://tesch1.gitlab.io/cppduals
|
||
//
|
||
// (c)2019 Michael Tesch. tesch1@gmail.com
|
||
//
|
||
// See https://gitlab.com/tesch1/cppduals/blob/master/LICENSE.txt for
|
||
// license information.
|
||
//
|
||
// This Source Code Form is subject to the terms of the Mozilla
|
||
// Public License v. 2.0. If a copy of the MPL was not distributed
|
||
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
||
|
||
#ifndef CPPDUALS_DUAL
|
||
#define CPPDUALS_DUAL
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
#define _DUALS_CONSTEXPR constexpr
|
||
|
||
///// FOR LIBCPP /////
|
||
#include <cstddef>
|
||
#if !defined(_LIBCPP_BEGIN_NAMESPACE_STD)
|
||
namespace std {
|
||
template <typename T> class complex;
|
||
}
|
||
#else // LIBCPP
|
||
_LIBCPP_BEGIN_NAMESPACE_STD
|
||
template <typename T> class complex;
|
||
_LIBCPP_END_NAMESPACE_STD
|
||
#endif // LIBCPP
|
||
|
||
#include <limits>
|
||
#include <type_traits>
|
||
|
||
#endif // PARSED_BY_DOXYGEN
|
||
|
||
#if !defined(CPPDUALS_IGNORE_COMPILER_VERSION) && !defined(_WIN32)
|
||
#if __cplusplus < 201103L
|
||
#error CPPDUALS needs at least a C++11 compliant compiler
|
||
#endif
|
||
#endif
|
||
|
||
/// Configure whether system has POSIX extern int signgam;
|
||
#ifndef CPPDUALS_HAVE_SIGNGAM
|
||
#ifndef _WIN32
|
||
#ifndef __CYGWIN__
|
||
#define CPPDUALS_HAVE_SIGNGAM 1
|
||
#endif
|
||
#endif
|
||
#endif
|
||
namespace duals {
|
||
|
||
/**
|
||
\file dual
|
||
\brief Dual number class.
|
||
|
||
\mainpage cppduals
|
||
|
||
\ref duals/dual is a single-header [Dual
|
||
number](https://en.wikipedia.org/wiki/Dual_number) C++ template
|
||
library that implements numbers of the type \f$(a + b \cdot
|
||
\epsilon)\f$ where \f$ \epsilon \ne 0\f$, and \f$\epsilon^2 = 0\f$.
|
||
|
||
`duals::dual<>` can be used for "automatic" differentiation, it can
|
||
also recursively nest with itself for higher orders of differentiation
|
||
(ie for the second derivative use `duals::dual<duals::dual<T>>`). It
|
||
can also be used to differentiate parts of complex functions as
|
||
`std::complex<duals::dual<T>>`. This file can be used stand-alone for
|
||
dual number support, but for Eigen vectorization support rather the
|
||
file \ref duals/dual_eigen should be included.
|
||
|
||
```
|
||
#include <duals/dual>
|
||
|
||
using namespace duals::literals;
|
||
|
||
template <class T> T f(T x) { return pow(x,pow(x,x)); }
|
||
template <class T> T df(T x) { return pow(x,-1. + x + pow(x,x)) * (1. + x*log(x) + x*pow(log(x),2.)); }
|
||
template <class T> T ddf(T x) { return (pow(x,pow(x,x)) * pow(pow(x,x - 1.) + pow(x,x)*log(x)*(log(x) + 1.), 2.) +
|
||
pow(x,pow(x,x)) * (pow(x,x - 1.) * log(x) +
|
||
pow(x,x - 1.) * (log(x) + 1.) +
|
||
pow(x,x - 1.) * ((x - 1.)/x + log(x)) +
|
||
pow(x,x) * log(x) * pow(log(x) + 1., 2.) )); }
|
||
|
||
int main()
|
||
{
|
||
std::cout << " f(2.) = " << f(2.) << "\n";
|
||
std::cout << " df(2.) = " << df(2.) << "\n";
|
||
std::cout << "ddf(2.) = " << ddf(2.) << "\n";
|
||
std::cout << " f(2+1_e) = " << f(2+1_e) << "\n";
|
||
std::cout << " f(2+1_e).dpart() = " << f(2+1_e).dpart() << "\n";
|
||
duals::hyperduald x(2+1_e,1+0_e);
|
||
std::cout << " f((2+1_e) + (1+0_e)_e).dpart().dpart() = " << f(x).dpart().dpart() << "\n";
|
||
}
|
||
```
|
||
Produces (notice the derivative in the dual-part):
|
||
```
|
||
f(2.) = 16
|
||
df(2.) = 107.11
|
||
ddf(2.) = 958.755
|
||
f(2+1_e) = (16+107.11_e)
|
||
f(2+1_e).dpart() = 107.11
|
||
f((2+1_e) + (1+0_e)_e).dpart().dpart() = 958.755
|
||
```
|
||
|
||
How this works can be seen by inspecting the infinite [Taylor
|
||
series](https://en.wikipedia.org/wiki/Taylor_series) expansion of a
|
||
function \f$ f(x) \f$ at \f$ a \f$ with \f$ x = a + b\epsilon \f$.
|
||
The series truncates itself due to the property \f$ \epsilon^2 = 0
|
||
\f$, leaving the function's derivative at \f$ a \f$ in the "dual part"
|
||
of the result (when \f$ b = 1 \f$):
|
||
|
||
\f[
|
||
\begin{split}
|
||
f(a + b \epsilon) &= f(a) + f'(a)(b \epsilon) + \frac{f''(a)}{2!}(b \epsilon)^2
|
||
+ \frac{f'''(a)}{3!}(b \epsilon)^3 + \ldots \\ \
|
||
&= f(a) + f'(a)(b \epsilon)
|
||
\end{split}
|
||
\f]
|
||
|
||
The class is contained in a single c++11 header `#include
|
||
<duals/dual>`, for extended [Eigen](http://eigen.tuxfamily.org)
|
||
support, instead include the header \ref duals/dual_eigen "#include <duals/dual_eigen>".
|
||
|
||
Type X in the templates below can be any value which can be
|
||
assigned to value_type.
|
||
|
||
Type X also indicates a limitation to dual numbers of the same depth
|
||
but (possibly) different value_type as `duals::dual<T>`. For example,
|
||
you can assign (or add/sub/mul/div) `duals::dual<float>` and
|
||
`duals::dual<double>`, but you can not assign
|
||
`duals::dual<duals::dual<float>>` to `duals::dual<float>`.
|
||
|
||
Here is a synopsis of the class:
|
||
|
||
```
|
||
namespace duals {
|
||
|
||
template<class T> class dual {
|
||
|
||
typedef T value_type;
|
||
|
||
dual(const & re = T(), const & du = T());
|
||
dual(const dual &);
|
||
template<class X> dual(const dual<X> &);
|
||
|
||
T rpart() const;
|
||
T dpart() const;
|
||
|
||
void rpart(T);
|
||
void dpart(T);
|
||
|
||
dual<T> operator-() const;
|
||
dual<T> operator+() const;
|
||
|
||
dual<T> & operator= (const T &);
|
||
dual<T> & operator+=(const T &);
|
||
dual<T> & operator-=(const T &);
|
||
dual<T> & operator*=(const T &);
|
||
dual<T> & operator/=(const T &);
|
||
|
||
dual<T> & operator=(const dual<T> &);
|
||
template<class X> dual<T> & operator= (const dual<X> &);
|
||
template<class X> dual<T> & operator+=(const dual<X> &);
|
||
template<class X> dual<T> & operator-=(const dual<X> &);
|
||
template<class X> dual<T> & operator*=(const dual<X> &);
|
||
template<class X> dual<T> & operator/=(const dual<X> &);
|
||
|
||
// The comparison operators are not strictly well-defined,
|
||
// they are implemented as comparison of the real part.
|
||
|
||
bool operator ==(const X &b) const;
|
||
bool operator !=(const X &b) const;
|
||
|
||
bool operator <(const X &b) const;
|
||
bool operator >(const X &b) const;
|
||
bool operator <=(const X &b) const;
|
||
bool operator >=(const X &b) const;
|
||
bool operator <(const dual<X> &b) const;
|
||
bool operator >(const dual<X> &b) const;
|
||
bool operator <=(const dual<X> &b) const;
|
||
bool operator >=(const dual<X> &b) const;
|
||
|
||
};
|
||
|
||
// Non-member functions:
|
||
|
||
T rpart(dual<T>) // Real part
|
||
T dpart(dual<T>) // Dual part
|
||
dual<T> dconj(dual<T>) // Dual-conjugate
|
||
|
||
// Transcendental functions
|
||
dual<T> exp(dual<T>)
|
||
dual<T> log(dual<T>)
|
||
dual<T> log10(dual<T>)
|
||
dual<T> log2(dual<T>)
|
||
dual<T> logb(dual<T> x)
|
||
dual<T> pow(dual<T>, U)
|
||
dual<T> pow(U, dual<T>)
|
||
dual<T> pow(dual<T>, dual<T>)
|
||
dual<T> sqrt(dual<T>)
|
||
dual<T> cbrt(dual<T>)
|
||
dual<T> sin(dual<T>)
|
||
dual<T> cos(dual<T>)
|
||
dual<T> tan(dual<T>)
|
||
dual<T> asin(dual<T>)
|
||
dual<T> acos(dual<T>)
|
||
dual<T> atan(dual<T>)
|
||
dual<T> atan2(dual<T>, dual<T>)
|
||
dual<T> hypot(dual<T>, dual<T>)
|
||
dual<T> scalbn(dual<T> arg, int ex)
|
||
// TODO:
|
||
dual<T> sinh(dual<T>)
|
||
dual<T> cosh(dual<T>)
|
||
dual<T> tanh(dual<T>)
|
||
dual<T> asinh(dual<T>)
|
||
dual<T> acosh(dual<T>)
|
||
dual<T> atanh(dual<T>)
|
||
|
||
dual<T> erf(dual<T>)
|
||
dual<T> erfc(dual<T>)
|
||
dual<T> tgamma(dual<T>)
|
||
dual<T> lgamma(dual<T>)
|
||
|
||
// Non-differentiable operations on the real part.
|
||
dual<T> abs(dual<T> d)
|
||
dual<T> fabs(dual<T> d)
|
||
dual<T> fmax(dual<T> x, dual<T> y)
|
||
dual<T> fmin(dual<T> x, dual<T> y)
|
||
dual<T> frexp(dual<T> arg, int* exp )
|
||
dual<T> ldexp(dual<T> arg, int exp )
|
||
dual<T> trunc(dual<T> d)
|
||
dual<T> floor(dual<T> d)
|
||
dual<T> ceil(dual<T> d)
|
||
dual<T> round(dual<T> d)
|
||
|
||
// floating point functions
|
||
int fpclassify(dual<T> d)
|
||
bool isfinite(dual<T> d)
|
||
bool isnormal(dual<T> d)
|
||
bool isinf(dual<T> d)
|
||
bool isnan(dual<T> d)
|
||
bool signbit(dual<T> d)
|
||
dual<T> copysign(dual<T> x, dual<T> y)
|
||
|
||
// Utility functions
|
||
dual<T> random(dual<T> low = {0,0}, dual<T> high = {1,0})
|
||
dual<T> random2(dual<T> low = {0,0}, dual<T> high = {1,0})
|
||
|
||
// Stream IO
|
||
template<T> operator>>(basic_istream<charT, traits> &, dual<T> &);
|
||
template<T> operator<<(basic_ostream<charT, traits> &, const dual<T> &);
|
||
|
||
}
|
||
```
|
||
|
||
Some useful typedefs:
|
||
|
||
```
|
||
typedef dual<float> dualf;
|
||
typedef dual<double> duald;
|
||
typedef dual<long double> dualld;
|
||
typedef dual<dualf> hyperdualf;
|
||
typedef dual<duald> hyperduald;
|
||
typedef dual<dualld> hyperdualld;
|
||
typedef std::complex<dualf> cdualf;
|
||
typedef std::complex<duald> cduald;
|
||
typedef std::complex<dualld> cdualld;
|
||
```
|
||
|
||
There are also literals for dual parts defined in the inline
|
||
namespace duals::literals. They are available with `using
|
||
namespace duals;` or `using namespace duals::literals`.
|
||
|
||
```
|
||
using namespace duals::literals;
|
||
dualf x = 3 + 4_ef;
|
||
duald y = 3 + 4_e;
|
||
dualld z = 3 + 4_el;
|
||
```
|
||
|
||
And in case you dislike iostreams, there are some formatters for the
|
||
[`{fmt}`](https://github.com/fmtlib/fmt) formatting library. These
|
||
are disabled by default, but can be enabled by `#define
|
||
CPPDUALS_LIBFMT` for the `dual<>` formatter, and/or `#define
|
||
CPPDUALS_LIBFMT_COMPLEX` for the std::complex<> formatter. There are
|
||
three custom formatting flags that control how the dual numbers are
|
||
printed, these must come before the normal `{fmt}` formatting spec
|
||
'$', '*', ','. For me these are about 3x faster than iostreams.
|
||
|
||
```
|
||
#define CPPDUALS_LIBFMT
|
||
#define CPPDUALS_LIBFMT_COMPLEX
|
||
#inlcude <duals/dual>
|
||
using namespace duals::literals;
|
||
...
|
||
string s = fmt::format("{:}", 1 + 2_e); // s = "(1.0+2.0_e)"
|
||
string s = fmt::format("{:g}", 1 + 2_e); // s = "(1+2_e)"
|
||
string s = fmt::format("{:*}", 1 + 2_e); // s = "(1.0+2.0*e)"
|
||
string s = fmt::format("{:,}", 1 + 2_e); // s = "(1.0,2.0)"
|
||
string s = fmt::format("{:*,g}", complexd(1,2_e)); // s = "((1)+(0,2)*i)"
|
||
```
|
||
|
||
The "archaic Greek epsilon" logo is from [Wikimedia
|
||
commons](https://commons.wikimedia.org/wiki/File:Greek_Epsilon_archaic.svg)
|
||
|
||
Some casual reading material:
|
||
|
||
- [ON QUATERNIONS, William Rowan Hamilton, Proceedings of the Royal Irish Academy, 3 (1847), pp. 1–16.](https://www.maths.tcd.ie/pub/HistMath/People/Hamilton/Quatern2/)
|
||
- [Basic Space-Time Transformations Expressed by Means of Two-Component Number Systems](https://doi.org/10.12693%2Faphyspola.86.291)
|
||
|
||
*/
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
template<class T> class dual;
|
||
#endif
|
||
|
||
/// Check if T is dual<>, match non-duals.
|
||
template <class T> struct is_dual : std::false_type {};
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
|
||
/// Check if T is dual<>, match dual<>.
|
||
template <class T> struct is_dual<dual<T> > : std::true_type {};
|
||
|
||
#endif
|
||
|
||
/// Check if T is std::complex<>, match non- std::complex<>.
|
||
template <class T> struct is_complex : std::false_type {};
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
|
||
/// Check if T is std::complex<>, match std::complex<>.
|
||
template <class T> struct is_complex<std::complex<T> > : std::true_type {};
|
||
|
||
#endif
|
||
|
||
/// Dual_traits helper class.
|
||
template <class T> struct dual_traits
|
||
{
|
||
/// Depth of T - for T=scalar this is 0. for dual_traits<double> it
|
||
/// is 1.
|
||
enum { depth = 0 }; // -Wenum-compare
|
||
|
||
/// The real storage type.
|
||
typedef T real_type;
|
||
};
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
|
||
/// dual_traits for dual<> types
|
||
template <class T> struct dual_traits<dual<T>>
|
||
{
|
||
/// Depth to which this dual<> type is nested. One (1) is a
|
||
/// first-level dual, whereas non-duals have a depth of 0.
|
||
enum { depth = dual_traits<T>::depth + 1 };
|
||
|
||
/// The real storage type.
|
||
typedef typename dual_traits<T>::real_type real_type;
|
||
};
|
||
|
||
template <class T> struct dual_traits<std::complex<dual<T>>>
|
||
{
|
||
/// complex<dual<T>> have the same 'depth' as their dual.
|
||
enum { depth = dual_traits<dual<T> >::depth };
|
||
|
||
/// The real storage type.
|
||
typedef typename dual_traits<T>::real_type real_type;
|
||
};
|
||
|
||
namespace detail {
|
||
|
||
template<class T>
|
||
struct Void { typedef void type; };
|
||
template<class T, class U = void>
|
||
struct has_member_type : std::false_type {};
|
||
template<class T>
|
||
struct has_member_type<T, typename Void<typename T::type>::type > : std::true_type {
|
||
struct wrap {
|
||
typedef typename T::type type;
|
||
typedef typename T::type ReturnType;
|
||
};
|
||
};
|
||
|
||
} // namespace detail
|
||
|
||
/// Promote two types - default according to common_type
|
||
template <class T, class U, class V = void> struct promote : std::common_type<T,U> {};
|
||
|
||
/// Can types A and B be promoted to a common type?
|
||
template<class A, class B> using can_promote = detail::has_member_type<promote<A,B>>;
|
||
|
||
// dual<T> dual<U>
|
||
// T == U
|
||
template <class T, class U>
|
||
struct promote<dual<T>, dual<U>,
|
||
typename std::enable_if<(can_promote<T,U>::value
|
||
&& (int)dual_traits<T>::depth == (int)dual_traits<U>::depth)>::type>
|
||
{
|
||
typedef dual<typename promote<U,T>::type> type;
|
||
};
|
||
// dual<T> U
|
||
// T >= U
|
||
// U is not complex
|
||
template <class T, class U>
|
||
struct promote<dual<T>, U,
|
||
typename std::enable_if<(can_promote<T,U>::value
|
||
&& (int)dual_traits<T>::depth >= (int)dual_traits<U>::depth
|
||
&& !is_complex<U>::value)>::type>
|
||
{
|
||
typedef dual<typename promote<U,T>::type> type;
|
||
};
|
||
// U dual<T>
|
||
// T >= U
|
||
// U is not complex
|
||
template <class T, class U>
|
||
struct promote<U, dual<T>,
|
||
typename std::enable_if<(can_promote<T,U>::value
|
||
&& (int)dual_traits<T>::depth >= (int)dual_traits<U>::depth
|
||
&& !is_complex<U>::value)>::type>
|
||
{
|
||
typedef dual<typename promote<U,T>::type> type;
|
||
};
|
||
// /////////////////////////////////////////////////
|
||
// complex<T> complex<U>
|
||
// T == U
|
||
template <class T, class U>
|
||
struct promote<std::complex<T>, std::complex<U>,
|
||
typename std::enable_if<(can_promote<T,U>::value
|
||
&& (is_dual<T>::value || is_dual<U>::value))>::type>
|
||
{
|
||
typedef std::complex<typename promote<T,U>::type> type;
|
||
};
|
||
// complex<T> U
|
||
// T >= U
|
||
// U is not complex
|
||
template <class T, class U>
|
||
struct promote<std::complex<T>, U,
|
||
typename std::enable_if<(can_promote<T,U>::value
|
||
&& (is_dual<T>::value || is_dual<U>::value)
|
||
&& !is_complex<U>::value)>::type>
|
||
{
|
||
typedef std::complex<typename promote<T,U>::type> type;
|
||
};
|
||
// U complex<T>
|
||
// T >= U
|
||
// U is not complex
|
||
template <class T, class U>
|
||
struct promote<U, std::complex<T>,
|
||
typename std::enable_if<(can_promote<T,U>::value
|
||
&& (is_dual<T>::value || is_dual<U>::value)
|
||
&& !is_complex<U>::value)>::type>
|
||
{
|
||
typedef std::complex<typename promote<T,U>::type> type;
|
||
};
|
||
// /////////////////////////////////////////////////
|
||
|
||
#endif // PARSED_BY_DOXYGEN
|
||
|
||
} // namespace duals
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
|
||
#define NOMACRO // thwart macroification
|
||
#ifdef EIGEN_PI
|
||
#define MY_PI EIGEN_PI
|
||
#else
|
||
#define MY_PI M_PI
|
||
#endif
|
||
|
||
#endif // PARSED_BY_DOXYGEN
|
||
|
||
#ifdef _LIBCPP_BEGIN_NAMESPACE_STD
|
||
_LIBCPP_BEGIN_NAMESPACE_STD
|
||
#else
|
||
namespace std {
|
||
#endif
|
||
|
||
#ifdef CPPDUALS_ENABLE_STD_IS_ARITHMETIC
|
||
|
||
/// Duals are as arithmetic as their value_type is arithmetic.
|
||
template <class T>
|
||
struct is_arithmetic<duals::dual<T>> : is_arithmetic<T> {};
|
||
|
||
#endif // CPPDUALS_ENABLE_IS_ARITHMETIC
|
||
|
||
/// Duals are compound types.
|
||
template <class T>
|
||
struct is_compound<duals::dual<T>> : true_type {};
|
||
|
||
// Modification of std::numeric_limits<> per
|
||
// C++03 17.4.3.1/1, and C++11 18.3.2.3/1.
|
||
template <class T>
|
||
struct numeric_limits<duals::dual<T> > : public numeric_limits<T> {
|
||
static constexpr bool is_specialized = true;
|
||
static constexpr duals::dual<T> min NOMACRO () { return numeric_limits<T>::min NOMACRO (); }
|
||
static constexpr duals::dual<T> lowest() { return numeric_limits<T>::lowest(); }
|
||
static constexpr duals::dual<T> max NOMACRO () { return numeric_limits<T>::max NOMACRO (); }
|
||
static constexpr duals::dual<T> epsilon() { return numeric_limits<T>::epsilon(); }
|
||
static constexpr duals::dual<T> round_error() { return numeric_limits<T>::round_error(); }
|
||
static constexpr duals::dual<T> infinity() { return numeric_limits<T>::infinity(); }
|
||
static constexpr duals::dual<T> quiet_NaN() { return numeric_limits<T>::quiet_NaN(); }
|
||
static constexpr duals::dual<T> signaling_NaN() { return numeric_limits<T>::signaling_NaN(); }
|
||
static constexpr duals::dual<T> denorm_min() { return numeric_limits<T>::denorm_min(); }
|
||
};
|
||
|
||
#ifdef _LIBCPP_BEGIN_NAMESPACE_STD
|
||
_LIBCPP_END_NAMESPACE_STD
|
||
#else
|
||
} // namespace std
|
||
#endif
|
||
|
||
namespace duals {
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
|
||
// T and X are wrapped in a dual<>
|
||
#define CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X) \
|
||
typename std::enable_if<(int)duals::dual_traits<X>::depth == \
|
||
(int)duals::dual_traits<T>::depth, int>::type = 0, \
|
||
typename std::enable_if<can_promote<T,X>::value,int>::type = 0
|
||
|
||
// Both T and U are wrapped in a dual<>
|
||
#define CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U) \
|
||
typename common_t = dual<typename duals::promote<T,U>::type>, \
|
||
typename std::enable_if<(int)duals::dual_traits<T>::depth == \
|
||
(int)duals::dual_traits<U>::depth, int>::type = 0, \
|
||
typename std::enable_if<can_promote<T,U>::value,int>::type = 0
|
||
|
||
// T is wrapped in a dual<>, U may or may not be.
|
||
#define CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U) \
|
||
typename common_t = typename duals::promote<dual<T>, U>::type, \
|
||
typename std::enable_if<((int)duals::dual_traits<T>::depth >= \
|
||
(int)duals::dual_traits<U>::depth), int>::type = 0, \
|
||
typename std::enable_if<can_promote<dual<T>,U>::value,int>::type = 0
|
||
|
||
// T is wrapped in complex<dual<>>
|
||
#define CPPDUALS_ENABLE_LEQ_DEPTH_AND_CX_COMMON_T(T,U) \
|
||
typename common_t = typename duals::promote<std::complex<dual<T> >,U>::type, \
|
||
typename std::enable_if<((int)duals::dual_traits<T>::depth >= \
|
||
(int)duals::dual_traits<U>::depth), int>::type = 0, \
|
||
typename std::enable_if<can_promote<std::complex<dual<T>>,U>::value,int>::type = 0
|
||
|
||
#define CPPDUALS_ENABLE_IF(...) typename std::enable_if< (__VA_ARGS__) , int>::type = 0
|
||
|
||
#endif
|
||
|
||
/*! \page user Background
|
||
TODO: Add text here...
|
||
*/
|
||
|
||
/// Abstract dual number class. Can nest with other dual numbers and
|
||
/// complex numbers.
|
||
template<class T>
|
||
class dual
|
||
{
|
||
public:
|
||
/// Class type of rpart() and dpart(). This type can be nested
|
||
/// dual<> or std::complex<>.
|
||
typedef T value_type;
|
||
|
||
private:
|
||
/// The real part.
|
||
value_type _real;
|
||
/// The dual part.
|
||
value_type _dual;
|
||
|
||
public:
|
||
|
||
/// Construct dual from optional real and dual parts.
|
||
constexpr
|
||
dual(const value_type re = value_type(), const value_type du = value_type())
|
||
: _real(re), _dual(du) {}
|
||
|
||
/// Copy construct from a dual of equal depth.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X),
|
||
CPPDUALS_ENABLE_IF(!is_complex<X>::value)>
|
||
dual(const dual<X> & x)
|
||
: _real(x.rpart()), _dual(x.dpart()) {}
|
||
|
||
/// Cast to a complex<dual<>> with real part equal to *this.
|
||
operator std::complex<dual<T>>() { return std::complex<dual<T>>(*this, (T)0); }
|
||
|
||
/// Explicit cast to an arithmetic type retains the rpart()
|
||
template <class X,
|
||
CPPDUALS_ENABLE_IF(std::is_arithmetic<X>::value && !is_dual<X>::value)>
|
||
explicit operator X() const { return X(_real); }
|
||
|
||
/// Get the real part.
|
||
T rpart() const { return _real; }
|
||
|
||
/// Get the dual part.
|
||
T dpart() const { return _dual; }
|
||
|
||
/// Set the real part.
|
||
void rpart(value_type re) { _real = re; }
|
||
|
||
/// Get the dual part.
|
||
void dpart(value_type du) { _dual = du; }
|
||
|
||
/// Unary negation
|
||
dual<T> operator-() const { return dual<T>(-_real, -_dual); }
|
||
|
||
/// Unary nothing
|
||
dual<T> operator+() const { return *this; }
|
||
|
||
/// Assignment of `value_type` assigns the real part and zeros the dual part.
|
||
dual<T> & operator= (const T & x) { _real = x; _dual = value_type(); return *this; }
|
||
|
||
/// Add a relatively-scalar to this dual.
|
||
dual<T> & operator+=(const T & x) { _real += x; return *this; }
|
||
|
||
/// Subtract a relatively-scalar from this dual.
|
||
dual<T> & operator-=(const T & x) { _real -= x; return *this; }
|
||
|
||
/// Multiply a relatively-scalar with this dual.
|
||
dual<T> & operator*=(const T & x) { _real *= x; _dual *= x; return *this; }
|
||
|
||
/// Divide this dual by relatively-scalar.
|
||
dual<T> & operator/=(const T & x) { _real /= x; _dual /= x; return *this; }
|
||
|
||
//dua & operator=(const dual & x) { _real = x.rpart(); _dual = x.dpart(); }
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
dual<T> & operator= (const dual<X> & x) { _real = x.rpart(); _dual = x.dpart(); return *this; }
|
||
|
||
/// Add a dual of the same depth to this dual.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
dual<T> & operator+=(const dual<X> & x) { _real += x.rpart(); _dual += x.dpart(); return *this; }
|
||
|
||
/// Subtract a dual of the same depth from this dual.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
dual<T> & operator-=(const dual<X> & x) { _real -= x.rpart(); _dual -= x.dpart(); return *this; }
|
||
|
||
/// Multiply this dual with a dual of same depth.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
dual<T> & operator*=(const dual<X> & x) {
|
||
_dual = _real * x.dpart() + _dual * x.rpart();
|
||
_real = _real * x.rpart();
|
||
return *this;
|
||
}
|
||
|
||
/// Divide this dual by another dual of the same or lower depth.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
dual<T> & operator/=(const dual<X> & x) {
|
||
_dual = (_dual * x.rpart() - _real * x.dpart()) / (x.rpart() * x.rpart());
|
||
_real = _real / x.rpart();
|
||
return *this;
|
||
}
|
||
|
||
//The following comparison operators are not strictly well-defined,
|
||
//they are implemented as comparison of the real parts.
|
||
#if 0
|
||
/// Compare this dual with another dual, comparing real parts for equality.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
bool operator ==(const dual<X> &b) const { return _real == b.rpart(); }
|
||
|
||
/// Compare this dual with another dual, comparing real parts for inequality.
|
||
template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
|
||
bool operator !=(const dual<X> &b) const { return _real != b.rpart(); }
|
||
|
||
/// Compare real part against real part of b
|
||
template<class X, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,X)>
|
||
bool operator <(const dual<X> &b) const { return _real < b.rpart(); }
|
||
|
||
/// Compare real part against real part of b
|
||
template<class X, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,X)>
|
||
bool operator >(const dual<X> &b) const { return _real > b.rpart(); }
|
||
|
||
/// Compare real part against real part of b
|
||
template<class X, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,X)>
|
||
bool operator <=(const dual<X> &b) const { return _real <= b.rpart(); }
|
||
|
||
/// Compare real part against real part of b
|
||
template<class X, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,X)>
|
||
bool operator >=(const dual<X> &b) const { return _real >= b.rpart(); }
|
||
#endif
|
||
};
|
||
|
||
/// Get the dual's real part.
|
||
template <class T> T rpart(const dual<T> & x) { return x.rpart(); }
|
||
|
||
/// Get the dual's dual part.
|
||
template <class T> T dpart(const dual<T> & x) { return x.dpart(); }
|
||
|
||
/// R-part of complex<dual<>> is non-dual complex<> (not to be confused with real())
|
||
template <class T> std::complex<T> rpart(const std::complex<dual<T>> & x)
|
||
{ return std::complex<T>(x.real().rpart(), x.imag().rpart()); }
|
||
|
||
/// Dual part of complex<dual<>> is complex<>
|
||
template <class T> std::complex<T> dpart(const std::complex<dual<T>> & x)
|
||
{ return std::complex<T>(x.real().dpart(), x.imag().dpart()); }
|
||
|
||
/// Get a non-dual's real part.
|
||
template <class T,
|
||
CPPDUALS_ENABLE_IF((std::is_arithmetic<T>::value &&
|
||
!std::is_compound<T>::value) || is_complex<T>::value)>
|
||
T rpart(const T & x) { return x; }
|
||
|
||
/// Get a non-dual's dual part := zero.
|
||
template <class T,
|
||
CPPDUALS_ENABLE_IF((std::is_arithmetic<T>::value &&
|
||
!std::is_compound<T>::value) || is_complex<T>::value)>
|
||
T dpart(const T & ) { return T(0); }
|
||
|
||
#ifndef PARSED_BY_DOXYGEN
|
||
|
||
/// Dual +-*/ ops with another entity
|
||
#define CPPDUALS_BINARY_OP(op) \
|
||
template<class T, class U, CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U)> \
|
||
common_t \
|
||
operator op(const dual<T> & z, const dual<U> & w) { \
|
||
common_t x(z); \
|
||
return x op##= w; \
|
||
} \
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U), \
|
||
CPPDUALS_ENABLE_IF(!std::is_same<U,std::complex<dual<T>>>::value)> \
|
||
common_t \
|
||
operator op(const dual<T> & z, const U & w) { \
|
||
common_t x(z); \
|
||
return x op##= w; \
|
||
} \
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U), \
|
||
CPPDUALS_ENABLE_IF(!std::is_same<U,std::complex<dual<T>>>::value)> \
|
||
common_t \
|
||
operator op(const U & z, const dual<T> & w) { \
|
||
common_t x(z); \
|
||
return x op##= w; \
|
||
} \
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_CX_COMMON_T(T,U), \
|
||
CPPDUALS_ENABLE_IF(!std::is_same<U,std::complex<dual<T>>>::value)> \
|
||
common_t \
|
||
operator op(const std::complex<dual<T>> & z, const U & w) { \
|
||
common_t x(z); \
|
||
return x op##= w; \
|
||
} \
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_CX_COMMON_T(T,U), \
|
||
CPPDUALS_ENABLE_IF(!std::is_same<U,std::complex<dual<T>>>::value)> \
|
||
common_t \
|
||
operator op(const U & z, const std::complex<dual<T>> & w) { \
|
||
common_t x(z); \
|
||
return x op##= w; \
|
||
}
|
||
|
||
CPPDUALS_BINARY_OP(+)
|
||
CPPDUALS_BINARY_OP(-)
|
||
CPPDUALS_BINARY_OP(*)
|
||
CPPDUALS_BINARY_OP(/)
|
||
|
||
/// Dual compared to a non-complex lower rank thing
|
||
#define CPPDUALS_LHS_COMPARISON(op) \
|
||
template<class T, class U, CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U)> \
|
||
bool \
|
||
operator op(const dual<T> & a, const dual<U> & b) { \
|
||
return a.rpart() op b.rpart(); \
|
||
} \
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U), \
|
||
CPPDUALS_ENABLE_IF(!is_complex<U>::value)> \
|
||
bool \
|
||
operator op(const U & a, const dual<T> & b) { \
|
||
return a op b.rpart(); \
|
||
} \
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U), \
|
||
CPPDUALS_ENABLE_IF(!is_complex<U>::value)> \
|
||
bool \
|
||
operator op(const dual<T> & a, const U & b) { \
|
||
return a.rpart() op b; \
|
||
}
|
||
|
||
CPPDUALS_LHS_COMPARISON(<)
|
||
CPPDUALS_LHS_COMPARISON(>)
|
||
CPPDUALS_LHS_COMPARISON(<=)
|
||
CPPDUALS_LHS_COMPARISON(>=)
|
||
CPPDUALS_LHS_COMPARISON(==)
|
||
CPPDUALS_LHS_COMPARISON(!=)
|
||
|
||
#endif // PARSED_BY_DOXYGEN
|
||
|
||
/// Complex Conjugate of a dual is just the dual.
|
||
template<class T> dual<T> conj(const dual<T> & x) { return x; }
|
||
|
||
/// Conjugate a thing that's not dual and not complex -- it has no
|
||
/// complex value so just return it. This is different from
|
||
/// std::conj() which promotes the T to a std::complex<T>.
|
||
template<class T, CPPDUALS_ENABLE_IF(!is_dual<T>::value &&
|
||
!is_complex<T>::value &&
|
||
std::is_arithmetic<T>::value)>
|
||
T conj(const T & x) { return x; }
|
||
|
||
/// Dual Conjugate, such that x*dconj(x) = rpart(x)^2. Different
|
||
/// function name from complex conjugate conj().
|
||
template<class T> dual<T> dconj(const dual<T> & x) {
|
||
return dual<T>(x.rpart(), - x.dpart());
|
||
}
|
||
|
||
/// Dual Conjugate of a complex, such that x*dconj(x) = rpart(x)^2.
|
||
/// Different function name from complex conjugate conj().
|
||
template<class T> std::complex<T> dconj(const std::complex<T> & x) {
|
||
return std::complex<T>(dconj(x.real()), dconj(x.imag()));
|
||
}
|
||
|
||
/// Conjugate a thing that's not dual and not complex.
|
||
template<class T, CPPDUALS_ENABLE_IF(!is_dual<T>::value &&
|
||
!is_complex<T>::value &&
|
||
std::is_arithmetic<T>::value)>
|
||
T dconj(const T & x) { return x; }
|
||
|
||
namespace detail {
|
||
template<class T> dual<T> exp_impl(const dual<T> & x);
|
||
template<class T> dual<T> log_impl(const dual<T> & x);
|
||
template<class T> dual<T> log10_impl(const dual<T> & x);
|
||
template<class T> dual<T> log2_impl(const dual<T> & x);
|
||
template<class T> dual<T> logb_impl(const dual<T> & x);
|
||
template<typename T, typename U, typename common_t>
|
||
common_t pow_dual_scalar(const dual<T>& x, const U& y);
|
||
template<typename T, typename U, typename common_t>
|
||
common_t pow_scalar_dual(const U& x, const dual<T>& y);
|
||
template<typename T>
|
||
std::complex<dual<T>> pow_complex_dual(const std::complex<dual<T>>& x, const dual<T>& y);
|
||
template<typename T>
|
||
std::complex<dual<T>> pow_dual_complex(const dual<T>& x, const std::complex<dual<T>>& y);
|
||
template<typename T, typename U, typename common_t>
|
||
common_t pow_dual_dual(const dual<T>& x, const dual<U>& y);
|
||
template<class T> dual<T> sqrt_impl(const dual<T> & x);
|
||
template<class T> dual<T> cbrt_impl(const dual<T> & x);
|
||
template<class T> dual<T> sin_impl(const dual<T> & x);
|
||
template<class T> dual<T> cos_impl(const dual<T> & x);
|
||
template<class T> dual<T> tan_impl(const dual<T> & x);
|
||
template<class T> dual<T> asin_impl(const dual<T> & x);
|
||
template<class T> dual<T> acos_impl(const dual<T> & x);
|
||
template<class T> dual<T> atan_impl(const dual<T> & x);
|
||
template<class T, class U, typename common_t>
|
||
common_t atan2_dual_scalar(const dual<T> & y, const U & x);
|
||
template<class T, class U, typename common_t>
|
||
common_t atan2_scalar_dual(const U & y, const dual<T> & x);
|
||
template<class T> dual<T> atan2_dual_dual(const dual<T> & y, const dual<T> & x);
|
||
template<class T> dual<T> hypot_impl(const dual<T> & x, const dual<T> & y);
|
||
template<class T> dual<T> scalbn_impl(const dual<T> & arg, int ex);
|
||
template<class T> dual<T> sinh_impl(const dual<T> & x);
|
||
template<class T> dual<T> cosh_impl(const dual<T> & x);
|
||
template<class T> dual<T> tanh_impl(const dual<T> & x);
|
||
template<class T> dual<T> asinh_impl(const dual<T> & x);
|
||
template<class T> dual<T> acosh_impl(const dual<T> & x);
|
||
template<class T> dual<T> atanh_impl(const dual<T> & x);
|
||
template<class T> dual<T> erf_impl(const dual<T> & x);
|
||
template<class T> dual<T> erfc_impl(const dual<T> & x);
|
||
template<class T> dual<T> tgamma_impl(const dual<T> & x);
|
||
template<class T> dual<T> lgamma_impl(const dual<T> & x);
|
||
template<class T> dual<T> abs_impl(const dual<T> & x);
|
||
template<class T> dual<T> fabs_impl(const dual<T> & x);
|
||
template<class T> dual<T> fmax_impl(const dual<T> & x, const dual<T> & y);
|
||
template<class T> dual<T> fmin_impl(const dual<T> & x, const dual<T> & y);
|
||
template<class T> dual<T> frexp_impl(const dual<T> & x, int* exp);
|
||
template<class T> dual<T> ldexp_impl(const dual<T> & x, int exp);
|
||
template<class T> dual<T> trunc_impl(const dual<T> & x);
|
||
template<class T> dual<T> floor_impl(const dual<T> & x);
|
||
template<class T> dual<T> ceil_impl(const dual<T> & x);
|
||
template<class T> dual<T> round_impl(const dual<T> & x);
|
||
template<class T> int fpclassify_impl(const dual<T> & x);
|
||
template<class T> bool isfinite_impl(const dual<T> & x);
|
||
template<class T> bool isnormal_impl(const dual<T> & x);
|
||
template<class T> bool isinf_impl(const dual<T> & x);
|
||
template<class T> bool isnan_impl(const dual<T> & x);
|
||
template<class T> bool signbit_impl(const dual<T> & x);
|
||
template<class T> dual<T> copysign_impl(const dual<T> & x, const dual<T> & y);
|
||
template<class T> dual<T> log1p_impl(const dual<T> & x);
|
||
template<class T> dual<T> expm1_impl(const dual<T> & x);
|
||
template<class T> dual<T> random_impl(const dual<T> & low, const dual<T> & high);
|
||
template<class T> dual<T> random2_impl(const dual<T> & low, const dual<T> & high);
|
||
}
|
||
|
||
/// /// /// /// /// /// /// ///
|
||
/// Prototypes
|
||
///
|
||
|
||
// Public interfaces that forward to implementations
|
||
template<class T> dual<T> exp(const dual<T> & x) { return detail::exp_impl(x); }
|
||
template<class T> dual<T> log(const dual<T> & x) { return detail::log_impl(x); }
|
||
template<class T> dual<T> log10(const dual<T> & x) { return detail::log10_impl(x); }
|
||
template<class T> dual<T> log2(const dual<T> & x) { return detail::log2_impl(x); }
|
||
template<class T> dual<T> logb(const dual<T> & x) { return detail::logb_impl(x); }
|
||
template<typename T, typename U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U), CPPDUALS_ENABLE_IF(!is_complex<U>::value)>
|
||
common_t pow(const dual<T>& x, const U& y) { return detail::pow_dual_scalar<T,U,common_t>(x, y); }
|
||
template<typename T, typename U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U), CPPDUALS_ENABLE_IF(!is_complex<U>::value)>
|
||
common_t pow(const U& x, const dual<T>& y) { return detail::pow_scalar_dual<T,U,common_t>(x, y); }
|
||
template <typename T>
|
||
std::complex<dual<T>> pow(const dual<T> & x, const std::complex<dual<T>> & y) { return detail::pow_dual_complex<T>(x, y); }
|
||
template <typename T>
|
||
std::complex<dual<T>> pow(const std::complex<dual<T>> & x, const dual<T>& y) { return detail::pow_complex_dual<T>(x, y); }
|
||
template <typename T, typename U, CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U)>
|
||
common_t pow(const dual<T>& x, const dual<U>& y) { return detail::pow_dual_dual<T,U,common_t>(x, y); }
|
||
template<class T> dual<T> sqrt(const dual<T> & x) { return detail::sqrt_impl(x); }
|
||
template<class T> dual<T> cbrt(const dual<T> & x) { return detail::cbrt_impl(x); }
|
||
template<class T> dual<T> sin(const dual<T> & x) { return detail::sin_impl(x); }
|
||
template<class T> dual<T> cos(const dual<T> & x) { return detail::cos_impl(x); }
|
||
template<class T> dual<T> tan(const dual<T> & x) { return detail::tan_impl(x); }
|
||
template<class T> dual<T> asin(const dual<T> & x) { return detail::asin_impl(x); }
|
||
template<class T> dual<T> acos(const dual<T> & x) { return detail::acos_impl(x); }
|
||
template<class T> dual<T> atan(const dual<T> & x) { return detail::atan_impl(x); }
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>
|
||
common_t atan2(const dual<T> & y, const U & x) { return detail::atan2_dual_scalar<T,U,common_t>(y, x); }
|
||
template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>
|
||
common_t atan2(const U & y, const dual<T> & x) { return detail::atan2_scalar_dual<T,U,common_t>(y, x); }
|
||
template<class T> dual<T> atan2(const dual<T> & y, const dual<T> & x) { return detail::atan2_dual_dual(y, x); }
|
||
template<class T> dual<T> hypot(const dual<T> & x, const dual<T> & y) { return detail::hypot_impl(x, y); }
|
||
template<class T> dual<T> scalbn(const dual<T> & arg, int ex) { return detail::scalbn_impl(arg, ex); }
|
||
template<class T> dual<T> sinh(const dual<T> & x) { return detail::sinh_impl(x); }
|
||
template<class T> dual<T> cosh(const dual<T> & x) { return detail::cosh_impl(x); }
|
||
template<class T> dual<T> tanh(const dual<T> & x) { return detail::tanh_impl(x); }
|
||
template<class T> dual<T> asinh(const dual<T> & x) { return detail::asinh_impl(x); }
|
||
template<class T> dual<T> acosh(const dual<T> & x) { return detail::acosh_impl(x); }
|
||
template<class T> dual<T> atanh(const dual<T> & x) { return detail::atanh_impl(x); }
|
||
/// The error function. Make sure to `#include <math.h>` before
|
||
/// `#include <duals/dual>` to use this function.
|
||
template<class T> dual<T> erf(const dual<T> & x) { return detail::erf_impl(x); }
|
||
template<class T> dual<T> erfc(const dual<T> & x) { return detail::erfc_impl(x); }
|
||
template<class T> dual<T> tgamma(const dual<T> & x) { return detail::tgamma_impl(x); }
|
||
template<class T> dual<T> lgamma(const dual<T> & x) { return detail::lgamma_impl(x); }
|
||
template<class T> dual<T> abs(const dual<T> & x) { return detail::abs_impl(x); }
|
||
template<class T> dual<T> fabs(const dual<T> & x) { return detail::fabs_impl(x); }
|
||
template<class T> dual<T> fmax(const dual<T> & x, const dual<T> & y) { return detail::fmax_impl(x, y); }
|
||
template<class T> dual<T> fmin(const dual<T> & x, const dual<T> & y) { return detail::fmin_impl(x, y); }
|
||
|
||
template<class T> dual<T> frexp(const dual<T> & x, int* exp) { return detail::frexp_impl(x, exp); }
|
||
template<class T> dual<T> ldexp(const dual<T> & x, int exp) { return detail::ldexp_impl(x, exp); }
|
||
//template<class T> dual<T> trunc(const dual<T> & x) { return detail::trunc_impl(x); }
|
||
template<class T> dual<T> trunc(const dual<T> & x) { return detail::trunc_impl(x); }
|
||
template<class T> dual<T> floor(const dual<T> & x) { return detail::floor_impl(x); }
|
||
template<class T> dual<T> ceil(const dual<T> & x) { return detail::ceil_impl(x); }
|
||
template<class T> dual<T> round(const dual<T> & x) { return detail::round_impl(x); }
|
||
|
||
template<class T> int fpclassify(const dual<T> & x) { return detail::fpclassify_impl(x); }
|
||
template<class T> bool isfinite(const dual<T> & x) { return detail::isfinite_impl(x); }
|
||
template<class T> bool isnormal(const dual<T> & x) { return detail::isnormal_impl(x); }
|
||
template<class T> bool isinf(const dual<T> & x) { return detail::isinf_impl(x); }
|
||
template<class T> bool isnan(const dual<T> & x) { return detail::isnan_impl(x); }
|
||
template<class T> bool signbit(const dual<T> & x) { return detail::signbit_impl(x); }
|
||
template<class T> dual<T> copysign(const dual<T> & x, const dual<T> & y) { return detail::copysign_impl(x, y); }
|
||
|
||
//template<class T> dual<T> random(const dual<T> & low, const dual<T> & high) { return detail::random_impl(low, high); }
|
||
//template<class T> dual<T> random2(const dual<T> & low, const dual<T> & high) { return detail::random2_impl(low, high); }
|
||
|
||
namespace utils {
|
||
template <typename T> int sgn(T val) {
|
||
return (T(0) < val) - (val < T(0));
|
||
}
|
||
}
|
||
|
||
template<class T> dual<T> log1p(const dual<T> & x) {return detail::log1p_impl(x);}
|
||
template<class T> dual<T> expm1(const dual<T> & x) {return detail::expm1_impl(x);}
|
||
|
||
#if __cpp_user_defined_literals >= 200809
|
||
/// Dual number literals in namespace duals::literals
|
||
inline namespace literals
|
||
{
|
||
using duals::dual;
|
||
constexpr dual<float> operator""_ef(long double du)
|
||
{
|
||
return { 0.0f, static_cast<float>(du) };
|
||
}
|
||
constexpr dual<float> operator""_ef(unsigned long long du)
|
||
{
|
||
return { 0.0f, static_cast<float>(du) };
|
||
}
|
||
constexpr dual<double> operator""_e(long double du)
|
||
{
|
||
return { 0.0, static_cast<double>(du) };
|
||
}
|
||
constexpr dual<double> operator""_e(unsigned long long du)
|
||
{
|
||
return { 0.0, static_cast<double>(du) };
|
||
}
|
||
constexpr dual<long double> operator""_el(long double du)
|
||
{
|
||
return { 0.0l, du };
|
||
}
|
||
constexpr dual<long double> operator""_el(unsigned long long du)
|
||
{
|
||
return { 0.0l, static_cast<long double>(du) };
|
||
}
|
||
} // namespace literals
|
||
#endif
|
||
|
||
// shortcut type names
|
||
typedef dual<float> dualf;
|
||
typedef dual<double> duald;
|
||
typedef dual<long double> dualld;
|
||
typedef dual<dualf> hyperdualf;
|
||
typedef dual<duald> hyperduald;
|
||
typedef dual<dualld> hyperdualld;
|
||
|
||
} // namespace duals
|
||
|
||
/// /// /// /// /// /// /// /// /// /// /// /// /// /// ///
|
||
/// Insert some stuff into std::
|
||
/// /// /// /// /// /// /// /// /// /// /// /// /// /// ///
|
||
|
||
#ifdef _LIBCPP_BEGIN_NAMESPACE_STD
|
||
_LIBCPP_BEGIN_NAMESPACE_STD
|
||
#else
|
||
namespace std {
|
||
#endif
|
||
|
||
// special needs LIBCPP
|
||
#if defined(_LIBCPP_BEGIN_NAMESPACE_STD)
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_copysign(duals::dual<DT> a, duals::dual<DT> b);
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_fabs(duals::dual<DT> a);
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isfinite(duals::dual<DT> a);
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isinf(duals::dual<DT> a);
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_fmax(duals::dual<DT> a, duals::dual<DT> b);
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_scalbn(duals::dual<DT> a, int exp);
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_logb(duals::dual<DT> a);
|
||
#endif
|
||
|
||
// export functions to the std namespace
|
||
using duals::exp;
|
||
using duals::log;
|
||
using duals::log10;
|
||
using duals::log2;
|
||
using duals::logb;
|
||
using duals::pow;
|
||
using duals::sqrt;
|
||
using duals::cbrt;
|
||
using duals::sin;
|
||
using duals::cos;
|
||
using duals::tan;
|
||
using duals::asin;
|
||
using duals::acos;
|
||
using duals::atan;
|
||
using duals::atan2;
|
||
using duals::hypot;
|
||
using duals::scalbn;
|
||
using duals::sinh;
|
||
using duals::cosh;
|
||
using duals::tanh;
|
||
using duals::asinh;
|
||
using duals::acosh;
|
||
using duals::atanh;
|
||
using duals::erf;
|
||
using duals::erfc;
|
||
using duals::tgamma;
|
||
using duals::lgamma;
|
||
using duals::abs;
|
||
using duals::fabs;
|
||
using duals::fmax;
|
||
using duals::fmin;
|
||
using duals::frexp;
|
||
using duals::ldexp;
|
||
using duals::trunc;
|
||
using duals::floor;
|
||
using duals::ceil;
|
||
using duals::round;
|
||
using duals::fpclassify;
|
||
using duals::isfinite;
|
||
using duals::isnormal;
|
||
using duals::isinf;
|
||
using duals::isnan;
|
||
using duals::signbit;
|
||
using duals::copysign;
|
||
//using duals::random;
|
||
//using duals::random2;
|
||
|
||
// TODO: figure out if still needed
|
||
namespace __math {
|
||
using duals::isinf;
|
||
}
|
||
|
||
#ifdef _LIBCPP_END_NAMESPACE_STD
|
||
_LIBCPP_END_NAMESPACE_STD
|
||
#else
|
||
} // namespace std
|
||
#endif
|
||
|
||
/// /// /// /// /// /// /// /// /// /// /// /// /// /// ///
|
||
/// Implementations
|
||
/// /// /// /// /// /// /// /// /// /// /// /// /// /// ///
|
||
|
||
#include <cmath>
|
||
#include <cerrno>
|
||
|
||
#ifdef _LIBCPP_BEGIN_NAMESPACE_STD
|
||
_LIBCPP_BEGIN_NAMESPACE_STD
|
||
#else
|
||
namespace std {
|
||
#endif
|
||
|
||
#if defined(_LIBCPP_BEGIN_NAMESPACE_STD)
|
||
// special needs LIBCPP
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_copysign(duals::dual<DT> a, duals::dual<DT> b) { return duals::copysign(a, b); }
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_fabs(duals::dual<DT> a) { return duals::fabs(a); }
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isfinite(duals::dual<DT> a) { return duals::isfinite(a); }
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isinf(duals::dual<DT> a) { return duals::isinf(a); }
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_fmax(duals::dual<DT> a, duals::dual<DT> b) { return duals::fmax(a, b); }
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_scalbn(duals::dual<DT> a, int exp) { return duals::scalbn(a, exp); }
|
||
template <class DT> _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual<DT> __constexpr_logb(duals::dual<DT> a) { return duals::logb(a); }
|
||
#endif
|
||
|
||
#ifdef _LIBCPP_END_NAMESPACE_STD
|
||
_LIBCPP_END_NAMESPACE_STD
|
||
#else
|
||
} // namespace std
|
||
#endif
|
||
|
||
namespace duals {
|
||
namespace detail {
|
||
|
||
/// Exponential e^x
|
||
template<class T> dual<T> exp_impl(const dual<T> & x) {
|
||
using std::exp;
|
||
T v = exp(x.rpart());
|
||
return dual<T>(v,
|
||
v * x.dpart());
|
||
}
|
||
|
||
/// Natural log ln(x)
|
||
template<class T> dual<T> log_impl(const dual<T> & x) {
|
||
using std::log;
|
||
T v = log(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else
|
||
return dual<T>(v, x.dpart() / x.rpart());
|
||
}
|
||
|
||
template<class T> dual<T> log10_impl(const dual<T> & x) {
|
||
using std::log;
|
||
return log(x) / log(static_cast<T>(10));
|
||
}
|
||
|
||
template<class T> dual<T> log2_impl(const dual<T> & x) {
|
||
using std::log;
|
||
return log(x) / log(static_cast<T>(2));
|
||
}
|
||
|
||
template<class T> duals::dual<T> logb_impl(const duals::dual<T> & x) {
|
||
using std::log;
|
||
using std::pow;
|
||
using std::abs;
|
||
using std::floor;
|
||
T v = floor(log(abs(x.rpart())) / log(std::numeric_limits<T>::radix));
|
||
return duals::dual<T>(v,
|
||
x.rpart() - pow(2,v) ? (T)0 : std::numeric_limits<T>::infinity());
|
||
}
|
||
|
||
// dual<T> ^ U implementation
|
||
template <typename T, typename U, typename common_t>
|
||
common_t pow_dual_scalar(const dual<T>& x, const U& y)
|
||
{
|
||
using std::pow;
|
||
typedef typename common_t::value_type V;
|
||
return common_t(
|
||
pow(x.rpart(), y),
|
||
x.dpart() == V(0) ? V(0) : x.dpart() * y * pow(x.rpart(), y - U(1)));
|
||
}
|
||
|
||
// U ^ dual<T> implementation
|
||
template <typename T, typename U, typename common_t>
|
||
common_t pow_scalar_dual(const U& x, const dual<T>& y)
|
||
{
|
||
return pow(common_t(x), y);
|
||
}
|
||
|
||
// U ^ dual<T> implementation
|
||
template <typename T>
|
||
std::complex<duals::dual<T>> pow_complex_dual(const std::complex<duals::dual<T>>& x, const dual<T>& y)
|
||
{
|
||
using dual_t = duals::dual<T>;
|
||
|
||
// 1. Compute the magnitude (r = |x|).
|
||
// r = sqrt( (Re x)^2 + (Im x)^2 ), where Re x and Im x are dual<T>.
|
||
dual_t rr = x.real() * x.real() + x.imag() * x.imag();
|
||
dual_t r = sqrt(rr); // Must have sqrt(dual_t) defined in your dual library
|
||
|
||
// 2. Compute the argument (theta = arg(x)).
|
||
// theta = atan2( Im x, Re x )
|
||
dual_t theta = atan2(x.imag(), x.real()); // dual-friendly atan2
|
||
|
||
// 3. Compute ln(x) = ln(r) + i * theta
|
||
// Then multiply by y: y * ln(x)
|
||
// realPart = y * ln(r)
|
||
// imagPart = y * theta
|
||
dual_t realPart = y * log(r);
|
||
dual_t imagPart = y * theta;
|
||
|
||
// 4. exponent = exp(realPart)
|
||
// Then real/imag of x^y = exponent * cos(imagPart), exponent *
|
||
// sin(imagPart)
|
||
dual_t exponent = exp(realPart);
|
||
|
||
dual_t realOut = exponent * cos(imagPart);
|
||
dual_t imagOut = exponent * sin(imagPart);
|
||
|
||
// 5. Return the result as a complex<dual<T>>.
|
||
return std::complex<dual_t>(realOut, imagOut);
|
||
}
|
||
|
||
template <typename T>
|
||
std::complex<duals::dual<T>> pow_dual_complex(
|
||
const dual<T>& x, const std::complex<duals::dual<T>>& y)
|
||
{
|
||
// For clarity, define short aliases:
|
||
using dual_t = duals::dual<T>;
|
||
|
||
// Extract the real and imaginary parts of y: both are dual<T>.
|
||
const dual_t& a = y.real(); // exponent's real part
|
||
const dual_t& b = y.imag(); // exponent's imaginary part
|
||
|
||
// 1. Compute ln(x). We assume x > 0 in the "real" sense.
|
||
// (In forward-mode dual, x can carry derivatives, but rpart(x) should be >
|
||
// 0
|
||
// or else you'd need branch handling for negative bases.)
|
||
dual_t ln_x = log(x);
|
||
|
||
// 2. Compute the real exponent = a * ln_x
|
||
// and the imaginary exponent = b * ln_x
|
||
dual_t realExp = a * ln_x; // a ln(x)
|
||
dual_t imagExp = b * ln_x; // b ln(x)
|
||
|
||
// 3. The magnitude is exp(a ln(x)):
|
||
dual_t magnitude = exp(realExp);
|
||
|
||
// 4. Then the final real and imaginary parts come from:
|
||
// Re(x^y) = magnitude * cos(b ln(x))
|
||
// Im(x^y) = magnitude * sin(b ln(x))
|
||
dual_t realPart = magnitude * cos(imagExp);
|
||
dual_t imagPart = magnitude * sin(imagExp);
|
||
|
||
// 5. Return the complex result.
|
||
return std::complex<dual_t>(realPart, imagPart);
|
||
}
|
||
|
||
// dual<T> ^ dual<U> implementation
|
||
template <typename T, typename U, typename common_t>
|
||
common_t pow_dual_dual(const dual<T>& f, const dual<U>& g)
|
||
{
|
||
using std::log;
|
||
using std::pow;
|
||
#if 1
|
||
using std::floor;
|
||
typedef typename common_t::value_type V;
|
||
common_t result;
|
||
|
||
if (f.rpart() == T(0) && g.rpart() >= U(1)) {
|
||
if (g.rpart() > U(1)) {
|
||
result = common_t(0);
|
||
}
|
||
else {
|
||
result = f;
|
||
}
|
||
}
|
||
else {
|
||
if (f.rpart() < T(0) && g.rpart() == floor(g.rpart())) {
|
||
V const tmp = g.rpart() * pow(f.rpart(), g.rpart() - U(1.0));
|
||
result = common_t(pow(f.rpart(), g.rpart()),
|
||
f.dpart() == T(0) ? T(0) : f.dpart() * tmp);
|
||
if (g.dpart() != U(0.0)) {
|
||
// Return a NaN when g.dpart() != 0.
|
||
result.dpart(std::numeric_limits<T>::quiet_NaN());
|
||
}
|
||
}
|
||
else {
|
||
// Handle the remaining cases. For cases 4,5,6,9 we allow the log()
|
||
// function to generate -HUGE_VAL or NaN, since those cases result in
|
||
// a nonfinite derivative.
|
||
V const tmp1 = pow(f.rpart(), g.rpart());
|
||
V const tmp2 = g.rpart() * pow(f.rpart(), g.rpart() - T(1.0));
|
||
V const tmp3 = tmp1 * log(f.rpart());
|
||
result = common_t(tmp1, f.dpart() == T(0) && g.dpart() == U(0)
|
||
? V(0)
|
||
: tmp2 * f.dpart() + tmp3 * g.dpart());
|
||
}
|
||
}
|
||
return result;
|
||
#else
|
||
T v = pow(f.rpart(), g.rpart());
|
||
return common_t(
|
||
v, pow(f.rpart(), g.rpart() - T(1)) *
|
||
(g.rpart() * f.dpart() + f.rpart() * log(f.rpart()) * g.dpart()));
|
||
#endif
|
||
}
|
||
|
||
template<class T> dual<T> sqrt_impl(const dual<T> & x) {
|
||
using std::sqrt;
|
||
T v = sqrt(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else
|
||
return dual<T>(v, x.dpart() / (T(2) * v) );
|
||
}
|
||
|
||
template<class T> dual<T> cbrt_impl(const dual<T> & x) {
|
||
using std::cbrt;
|
||
T v = cbrt(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else
|
||
return dual<T>(v, x.dpart() / (T(3) * v * v) );
|
||
}
|
||
|
||
template<class T> dual<T> sin_impl(const dual<T> & x) {
|
||
using std::sin;
|
||
using std::cos;
|
||
return dual<T>(sin(x.rpart()),
|
||
x.dpart() * cos(x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> cos_impl(const dual<T> & x) {
|
||
using std::cos;
|
||
using std::sin;
|
||
return dual<T>(cos(x.rpart()),
|
||
-sin(x.rpart()) * x.dpart());
|
||
}
|
||
|
||
template<class T> dual<T> tan_impl(const dual<T> & x) {
|
||
using std::tan;
|
||
T v = tan(x.rpart());
|
||
return dual<T>(v, x.dpart() * (v*v + 1));
|
||
}
|
||
|
||
template<class T> dual<T> asin_impl(const dual<T> & x) {
|
||
using std::asin;
|
||
using std::sqrt;
|
||
T v = asin(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else
|
||
return dual<T>(v, x.dpart() / sqrt(1 - x.rpart()*x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> acos_impl(const dual<T> & x) {
|
||
using std::acos;
|
||
using std::sqrt;
|
||
T v = acos(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else
|
||
return dual<T>(v, -x.dpart() / sqrt(1 - x.rpart()*x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> atan_impl(const dual<T> & x) {
|
||
using std::atan;
|
||
T v = atan(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else
|
||
return dual<T>(v, x.dpart() / (1 + x.rpart() * x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> atan2_dual_dual(const dual<T> & y, const dual<T> & x) {
|
||
using std::atan2;
|
||
T v = atan2(y.rpart(), x.rpart());
|
||
return dual<T>(v,
|
||
(x.rpart() * y.dpart() - y.rpart() * x.dpart()) /
|
||
(y.rpart() * y.rpart() + x.rpart() * x.rpart()));
|
||
}
|
||
|
||
template<class T, class U, typename common_t>
|
||
common_t
|
||
atan2_dual_scalar(const dual<T> & y, const U & x) {
|
||
using std::atan2;
|
||
T v = atan2(y.rpart(), x);
|
||
return dual<T>(v,
|
||
(x * y.dpart()) /
|
||
(y.rpart() * y.rpart() + x * x));
|
||
}
|
||
|
||
template<class T, class U, typename common_t>
|
||
common_t
|
||
atan2_scalar_dual(const U & y, const dual<T> & x) {
|
||
using std::atan2;
|
||
T v = atan2(y, x.rpart());
|
||
return dual<T>(v,
|
||
(-y * x.dpart()) /
|
||
(y * y + x.rpart() * x.rpart()));
|
||
}
|
||
|
||
template<class T> duals::dual<T> hypot_impl(const duals::dual<T> & x, const duals::dual<T> & y) {
|
||
using std::hypot;
|
||
return dual<T>(hypot(x.rpart(), y.rpart()),
|
||
(x.rpart() * y.dpart() + y.rpart() * x.dpart()) / hypot(x.rpart(), y.rpart()));
|
||
}
|
||
|
||
template<class T> duals::dual<T> scalbn_impl(const duals::dual<T> & arg, int ex) {
|
||
return arg * std::pow((T)2, ex);
|
||
}
|
||
|
||
template<class T> dual<T> sinh_impl(const dual<T> & x) {
|
||
using std::sinh;
|
||
using std::cosh;
|
||
return dual<T>(sinh(x.rpart()),
|
||
x.dpart() * cosh(x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> cosh_impl(const dual<T> & x) {
|
||
using std::cosh;
|
||
using std::sinh;
|
||
return dual<T>(cosh(x.rpart()),
|
||
x.dpart() * sinh(x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> tanh_impl(const dual<T> & x) {
|
||
using std::tanh;
|
||
return dual<T>(tanh(x.rpart()),
|
||
x.dpart() * (1 - tanh(x.rpart()) * tanh(x.rpart())));
|
||
}
|
||
|
||
template<class T> dual<T> asinh_impl(const dual<T> & x) {
|
||
using std::asinh;
|
||
using std::sqrt;
|
||
return dual<T>(asinh(x.rpart()),
|
||
x.dpart() / sqrt(1 + x.rpart() * x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> acosh_impl(const dual<T> & x) {
|
||
using std::acosh;
|
||
using std::sqrt;
|
||
return dual<T>(acosh(x.rpart()),
|
||
x.dpart() / sqrt(x.rpart() * x.rpart() - 1));
|
||
}
|
||
|
||
template<class T> dual<T> atanh_impl(const dual<T> & x) {
|
||
using std::atanh;
|
||
return dual<T>(atanh(x.rpart()),
|
||
x.dpart() / (1 - x.rpart() * x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> log1p_impl(const dual<T> & x) {
|
||
using std::log1p;
|
||
return dual<T>(log1p(x.rpart()),
|
||
x.dpart() / (1 + x.rpart()));
|
||
}
|
||
|
||
template<class T> dual<T> expm1_impl(const dual<T> & x) {
|
||
using std::expm1;
|
||
using std::exp;
|
||
return dual<T>(expm1(x.rpart()),
|
||
x.dpart() * exp(x.rpart()));
|
||
}
|
||
|
||
/// The error function. Make sure to `#include <math.h>` before
|
||
/// `#include <duals/dual>` to use this function.
|
||
template<class T> dual<T> erf_impl(const dual<T> & x) {
|
||
using std::erf;
|
||
using std::sqrt;
|
||
using std::pow;
|
||
using std::exp;
|
||
return dual<T>(erf(x.rpart()),
|
||
x.dpart() * T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))));
|
||
}
|
||
|
||
/// Error function complement (1 - erf()).
|
||
template<class T> dual<T> erfc_impl(const dual<T> & x) {
|
||
using std::erfc;
|
||
using std::sqrt;
|
||
using std::pow;
|
||
using std::exp;
|
||
return dual<T>(erfc(x.rpart()),
|
||
x.dpart() * -T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))));
|
||
}
|
||
|
||
/// Gamma function. Approximation of the dual part.
|
||
// TODO specialize for integers
|
||
template<class T> dual<T> tgamma_impl(const dual<T> & x) {
|
||
using std::tgamma;
|
||
T v = tgamma(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else {
|
||
|
||
int errno_saved = errno;
|
||
T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));
|
||
T w((tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h));
|
||
errno = errno_saved;
|
||
return dual<T>(v, x.dpart() * w);
|
||
}
|
||
}
|
||
|
||
/// Log of absolute value of gamma function. Approximation of the dual part.
|
||
template<class T> dual<T> lgamma_impl(const dual<T> & x) {
|
||
using std::lgamma;
|
||
T v = lgamma(x.rpart());
|
||
if (x.dpart() == T(0))
|
||
return v;
|
||
else {
|
||
#if CPPDUALS_HAVE_SIGNGAM
|
||
int signgam_saved = signgam;
|
||
#endif
|
||
int errno_saved = errno;
|
||
T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));
|
||
T w((lgamma(x.rpart()+h) - lgamma(x.rpart()-h))/(2*h));
|
||
#if CPPDUALS_HAVE_SIGNGAM
|
||
signgam = signgam_saved;
|
||
#endif
|
||
errno = errno_saved;
|
||
return dual<T>(v, x.dpart() * w);
|
||
}
|
||
}
|
||
template<class T> dual<T> abs_impl(const dual<T> & x) {
|
||
using std::abs;
|
||
return dual<T>(abs(x.rpart()),
|
||
x.dpart() * utils::sgn(x.rpart()));
|
||
}
|
||
template<class T> dual<T> fabs_impl(const dual<T> & x) {
|
||
using std::fabs;
|
||
return dual<T>(fabs(x.rpart()),
|
||
x.dpart() * utils::sgn(x.rpart()));
|
||
}
|
||
template<class T> duals::dual<T> fmax_impl(const duals::dual<T> & x, const duals::dual<T> & y) {
|
||
return x.rpart() > y.rpart() ? x : y;
|
||
}
|
||
template<class T> duals::dual<T> fmin_impl(const duals::dual<T> & x, const duals::dual<T> & y) {
|
||
return x.rpart() <= y.rpart() ? x : y;
|
||
}
|
||
template<class T> duals::dual<T> frexp_impl(const duals::dual<T> & arg, int* exp ) { using std::frexp; return frexp(arg.rpart(), exp); }
|
||
template<class T> duals::dual<T> ldexp_impl(const duals::dual<T> & arg, int exp ) { return arg * std::pow((T)2,exp); }
|
||
template<class T> duals::dual<T> trunc_impl(const duals::dual<T> & d) { using std::trunc; return trunc(d.rpart()); }
|
||
template<class T> dual<T> floor_impl(const dual<T> & x) {
|
||
using std::floor;
|
||
return dual<T>(floor(x.rpart()),
|
||
x.dpart() * utils::sgn(x.rpart()-floor(x.rpart())));
|
||
}
|
||
template<class T> duals::dual<T> ceil_impl(const duals::dual<T> & d) { using std::ceil; return ceil(d.rpart()); }
|
||
template<class T> duals::dual<T> round_impl(const duals::dual<T> & d) { using std::round; return round(d.rpart()); }
|
||
template<class T> int fpclassify_impl(const duals::dual<T> & d) { using std::fpclassify; return fpclassify(d.rpart()); }
|
||
template<class T> bool isfinite_impl(const duals::dual<T> & d) { using std::isfinite; return isfinite(d.rpart()); }
|
||
template<class T> bool isnormal_impl(const duals::dual<T> & d) { using std::isnormal; return isnormal(d.rpart()); }
|
||
template<class T> bool isinf_impl(const duals::dual<T> & d) { using std::isinf; return isinf(d.rpart()); }
|
||
template<class T> bool isnan_impl(const duals::dual<T> & d) { using std::isnan; return isnan(d.rpart()); }
|
||
template<class T> bool signbit_impl(const duals::dual<T> & d) { using std::signbit; return signbit(d.rpart()); }
|
||
template<class T> duals::dual<T> copysign_impl(const duals::dual<T> & x, const duals::dual<T> & y) {
|
||
using std::copysign;
|
||
T r = copysign(x.rpart(), y.rpart());
|
||
return duals::dual<T>(r, r == x.rpart() ? x.dpart() : -x.dpart());
|
||
}
|
||
|
||
} // namespace detail
|
||
} // namespace duals
|
||
|
||
#include <random>
|
||
namespace duals {
|
||
namespace randos {
|
||
|
||
// Random real value between a and b.
|
||
template<class T,
|
||
typename dist = std::uniform_real_distribution<T>,
|
||
CPPDUALS_ENABLE_IF(!is_complex<T>::value && !is_dual<T>::value)>
|
||
T random(T a = T(0), T b = T(1)) {
|
||
static std::default_random_engine generator;
|
||
dist distribution(a, b);
|
||
return distribution(generator);
|
||
}
|
||
|
||
template<class T,
|
||
typename dist = std::uniform_real_distribution<T>,
|
||
CPPDUALS_ENABLE_IF(!is_complex<T>::value && !is_dual<T>::value)>
|
||
T random2(T a = T(0), T b = T(1)) { return random<T,dist>(a,b); }
|
||
|
||
// Helper class for testing - also random value in dual part.
|
||
template <class DT,
|
||
CPPDUALS_ENABLE_IF(is_dual<DT>::value)>
|
||
DT random2(DT a = DT(0,0), DT b = DT(1,1)) {
|
||
using randos::random;
|
||
return DT(a.rpart() + random2<typename DT::value_type>() * (b.rpart() - a.rpart()),
|
||
a.dpart() + random2<typename DT::value_type>() * (b.dpart() - a.dpart()));
|
||
}
|
||
|
||
// Helper class for testing - also random value in dual part of the complex.
|
||
template<class CT,
|
||
CPPDUALS_ENABLE_IF(is_complex<CT>::value)>
|
||
CT random2(CT a = CT(0,0), CT b = CT(1,1)) {
|
||
return CT(a.real() + random2<typename CT::value_type>() * (b.real() - a.real()),
|
||
a.imag() + random2<typename CT::value_type>() * (b.imag() - a.imag()));
|
||
}
|
||
|
||
} // randos
|
||
|
||
/// Random real and dual parts, used by Eigen's Random(), by default
|
||
// the returned value has zero dual part and is in the range [0+0_e,
|
||
// 1+0_e].
|
||
template <class DT,
|
||
typename dist = std::uniform_real_distribution<typename DT::value_type>,
|
||
CPPDUALS_ENABLE_IF(is_dual<DT>::value)>
|
||
DT random(DT a = DT(0,0), DT b = DT(1,0)) {
|
||
using randos::random;
|
||
return DT(random<typename DT::value_type,dist>(a.rpart(),b.rpart()),
|
||
random<typename DT::value_type,dist>(a.dpart(),b.dpart()));
|
||
}
|
||
|
||
} // namespace duals
|
||
|
||
#include <complex>
|
||
|
||
namespace duals {
|
||
typedef std::complex<dualf> cdualf;
|
||
typedef std::complex<duald> cduald;
|
||
typedef std::complex<dualld> cdualld;
|
||
}
|
||
|
||
#include <iostream>
|
||
namespace duals {
|
||
|
||
/// Putto operator
|
||
template<class T, class _CharT, class _Traits>
|
||
std::basic_ostream<_CharT, _Traits> &
|
||
operator<<(std::basic_ostream<_CharT, _Traits> & os, const dual<T> & x)
|
||
{
|
||
std::basic_ostringstream<_CharT, _Traits> s;
|
||
s.flags(os.flags());
|
||
s.imbue(os.getloc());
|
||
s.precision(os.precision());
|
||
s << '(' << x.rpart()
|
||
<< (x.dpart() < 0 ? "" : "+")
|
||
<< x.dpart()
|
||
<< "_e" << (std::is_same<typename std::decay<T>::type,float>::value ? "f" :
|
||
std::is_same<typename std::decay<T>::type,double>::value ? "" :
|
||
std::is_same<typename std::decay<T>::type,long double>::value ? "l" : "")
|
||
<< ")";
|
||
return os << s.str();
|
||
}
|
||
|
||
/// Stream reader
|
||
template<class T, class CharT, class Traits>
|
||
std::basic_istream<CharT, Traits> &
|
||
operator>>(std::basic_istream<CharT, Traits> & is, dual<T> & x)
|
||
{
|
||
if (is.good()) {
|
||
ws(is);
|
||
if (is.peek() == CharT('(')) {
|
||
is.get();
|
||
T r;
|
||
is >> r;
|
||
if (!is.fail()) {
|
||
CharT c = is.peek();
|
||
if (c != CharT('_')) {
|
||
ws(is);
|
||
c = is.peek();
|
||
}
|
||
if (c == CharT('+') || c == CharT('-') || c == CharT('_')) {
|
||
if (c == CharT('+'))
|
||
is.get();
|
||
T d;
|
||
if (c == CharT('_')) {
|
||
d = r;
|
||
r = 0;
|
||
}
|
||
else
|
||
is >> d;
|
||
if (!is.fail()) {
|
||
ws(is);
|
||
c = is.peek();
|
||
if (c == CharT('_')) {
|
||
is.get();
|
||
c = is.peek();
|
||
if (c == CharT('e')) {
|
||
is.get();
|
||
c = is.peek();
|
||
if ((c == 'f' && !std::is_same<typename std::decay<T>::type,float>::value) ||
|
||
(c == 'l' && !std::is_same<typename std::decay<T>::type,long double>::value))
|
||
is.setstate(std::ios_base::failbit);
|
||
else {
|
||
if (c == 'f' || c == 'l')
|
||
is.get();
|
||
ws(is);
|
||
c = is.peek();
|
||
if (c == ')') {
|
||
is.get();
|
||
x = dual<T>(r, d);
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
else if (c == CharT(')')) {
|
||
is.get();
|
||
x = dual<T>(r, T(0));
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
else {
|
||
T r;
|
||
is >> r;
|
||
if (!is.fail())
|
||
x = dual<T>(r, T(0));
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
}
|
||
}
|
||
else
|
||
is.setstate(std::ios_base::failbit);
|
||
return is;
|
||
}
|
||
|
||
} // namespace duals
|
||
|
||
#ifdef CPPDUALS_LIBFMT
|
||
#include <fmt/format.h>
|
||
/// duals::dual<> Formatter for libfmt https://github.com/fmtlib/fmt
|
||
///
|
||
/// Formats a dual number (r,d) as (r+d_e), offering the same
|
||
/// formatting options as the underlying type - with the addition of
|
||
/// three optional format options, only one of which may appear
|
||
/// directly after the ':' in the format spec: '$', '*', and ','. The
|
||
/// '*' flag changes the separating _ to a *, producing (r+d*e), where
|
||
/// r and d are the formatted value_type values. The ',' flag simply
|
||
/// prints the real and dual parts separated by a comma. As a
|
||
/// concrete exmple, this formatter can produce either (3+5.4_e) or
|
||
/// (3+5.4*e) or (3,5.4) for a dual<double> using the specs {:g},
|
||
/// {:*g}, or {:,g}, respectively. When the '*' is NOT specified, the
|
||
/// output should be compatible with the input operator>> and the dual
|
||
/// literals below. (this implementation is a bit hacky - glad for
|
||
/// cleanups).
|
||
template <typename T, typename Char>
|
||
struct fmt::formatter<duals::dual<T>,Char> : public fmt::formatter<T,Char>
|
||
{
|
||
typedef fmt::formatter<T,Char> base;
|
||
enum style { expr, star, pair } style_ = expr;
|
||
detail::dynamic_format_specs<Char> specs_;
|
||
|
||
FMT_CONSTEXPR auto parse(parse_context<Char>& ctx) -> decltype(ctx.begin()) {
|
||
auto it = ctx.begin();
|
||
if (it != ctx.end()) {
|
||
switch (*it) {
|
||
case '$': style_ = style::expr; ctx.advance_to(++it); break;
|
||
case '*': style_ = style::star; ctx.advance_to(++it); break;
|
||
case ',': style_ = style::pair; ctx.advance_to(++it); break;
|
||
default: break;
|
||
}
|
||
}
|
||
return parse_format_specs(it, ctx.end(), specs_, ctx,
|
||
detail::type_constant<T, Char>::value);
|
||
}
|
||
|
||
template <typename FormatCtx>
|
||
auto format(const duals::dual<T> & x, FormatCtx & ctx) const -> decltype(ctx.out()) {
|
||
auto specs = specs_;
|
||
if (specs.dynamic()) {
|
||
detail::handle_dynamic_spec(specs.dynamic_width(), specs.width,
|
||
specs.width_ref, ctx);
|
||
detail::handle_dynamic_spec(specs.dynamic_precision(), specs.precision,
|
||
specs.precision_ref, ctx);
|
||
}
|
||
|
||
auto out = ctx.out();
|
||
*out++ = '(';
|
||
if (style_ == style::pair) {
|
||
out = detail::write<Char>(out, x.rpart(), specs, ctx.locale());
|
||
*out++ = ',';
|
||
out = detail::write<Char>(out, x.dpart(), specs, ctx.locale());
|
||
} else {
|
||
if (x.rpart() || !x.dpart())
|
||
out = detail::write<Char>(out, x.rpart(), specs, ctx.locale());
|
||
if (x.dpart()) {
|
||
if (x.rpart() && x.dpart() >= 0)
|
||
specs.set_sign(sign::plus);
|
||
out = detail::write<Char>(out, x.dpart(), specs, ctx.locale());
|
||
if (style_ == style::star) {
|
||
*out++ = '*';
|
||
*out++ = 'e';
|
||
} else {
|
||
*out++ = '_';
|
||
*out++ = 'e';
|
||
}
|
||
if (std::is_same<typename std::decay<T>::type,float>::value) *out++ = 'f';
|
||
if (std::is_same<typename std::decay<T>::type,long double>::value) *out++ = 'l';
|
||
}
|
||
}
|
||
*out++ = ')';
|
||
return out;
|
||
}
|
||
};
|
||
#endif // CPPDUALS_LIBFMT
|
||
|
||
#ifdef CPPDUALS_LIBFMT_COMPLEX
|
||
#ifndef CPPDUALS_LIBFMT
|
||
#include <fmt/format.h>
|
||
#endif
|
||
/// std::complex<> Formatter for libfmt https://github.com/fmtlib/fmt
|
||
///
|
||
/// libfmt does not provide a formatter for std::complex<>, although
|
||
/// one is proposed for c++20. Anyway, at the expense of a k or two,
|
||
/// you can define CPPDUALS_LIBFMT_COMPLEX and get this one.
|
||
///
|
||
/// The standard iostreams formatting of complex numbers is (a,b),
|
||
/// where a and b are the real and imaginary parts. This formats a
|
||
/// complex number (a+bi) as (a+bi), offering the same formatting
|
||
/// options as the underlying type - with the addition of three
|
||
/// optional format options, only one of which may appear directly
|
||
/// after the ':' in the format spec (before any fill or align): '$'
|
||
/// (the default if no flag is specified), '*', and ','. The '*' flag
|
||
/// adds a * before the 'i', producing (a+b*i), where a and b are the
|
||
/// formatted value_type values. The ',' flag simply prints the real
|
||
/// and complex parts separated by a comma (same as iostreams' format).
|
||
/// As a concrete exmple, this formatter can produce either (3+5.4i)
|
||
/// or (3+5.4*i) or (3,5.4) for a complex<double> using the specs {:g}
|
||
/// | {:$g}, {:*g}, or {:,g}, respectively. (this implementation is a
|
||
/// bit hacky - glad for cleanups).
|
||
///
|
||
template <typename T, typename Char>
|
||
struct fmt::formatter<std::complex<T>,Char> : public fmt::formatter<T,Char>
|
||
{
|
||
typedef fmt::formatter<T,Char> base;
|
||
enum style { expr, star, pair } style_ = expr;
|
||
detail::dynamic_format_specs<Char> specs_;
|
||
|
||
FMT_CONSTEXPR auto parse(parse_context<Char>& ctx) -> decltype(ctx.begin()) {
|
||
auto it = ctx.begin();
|
||
if (it != ctx.end()) {
|
||
switch (*it) {
|
||
case '$': style_ = style::expr; ctx.advance_to(++it); break;
|
||
case '*': style_ = style::star; ctx.advance_to(++it); break;
|
||
case ',': style_ = style::pair; ctx.advance_to(++it); break;
|
||
default: break;
|
||
}
|
||
}
|
||
// Let base parser handle any nested format chars first
|
||
auto end = base::parse(ctx);
|
||
// Now parse the format specs for this formatter
|
||
parse_format_specs(end, ctx.end(), specs_, ctx,
|
||
detail::type_constant<T, Char>::value);
|
||
return end;
|
||
}
|
||
|
||
template <typename FormatCtx>
|
||
auto format(const std::complex<T> & x, FormatCtx & ctx) const -> decltype(ctx.out()) {
|
||
format_to(ctx.out(), "(");
|
||
if (style_ == style::pair) {
|
||
base::format(x.real(), ctx);
|
||
format_to(ctx.out(), ",");
|
||
base::format(x.imag(), ctx);
|
||
} else {
|
||
if (x.real() || !x.imag())
|
||
base::format(x.real(), ctx);
|
||
if (x.imag()) {
|
||
// The base formatter will handle the negative sign
|
||
if (x.real() && x.imag() >= 0 && sign::plus != specs_.sign()) {
|
||
format_to(ctx.out(), "+");
|
||
}
|
||
base::format(x.imag(), ctx);
|
||
if (style_ == style::star)
|
||
format_to(ctx.out(), "*i");
|
||
else
|
||
format_to(ctx.out(), "i");
|
||
if (std::is_same<typename std::decay<T>::type,float>::value) format_to(ctx.out(), "f");
|
||
if (std::is_same<typename std::decay<T>::type,long double>::value) format_to(ctx.out(), "l");
|
||
}
|
||
}
|
||
return format_to(ctx.out(), ")");
|
||
}
|
||
};
|
||
#endif // CPPDUALS_LIBFMT_COMPLEX
|
||
|
||
#endif // CPPDUALS_DUAL
|