diff --git a/src/include/cppduals/.appveyor.yml b/src/include/cppduals/.appveyor.yml index 1bc8991f3..54587d2fd 100644 --- a/src/include/cppduals/.appveyor.yml +++ b/src/include/cppduals/.appveyor.yml @@ -3,9 +3,8 @@ clone_folder: c:\projects\cppduals clone_depth: 3 image: -#- Visual Studio 2013 -#- Visual Studio 2015 - Visual Studio 2017 +#- Visual Studio 2022 configuration: - Release @@ -14,12 +13,10 @@ configuration: # Do not build feature branch with open Pull Requests skip_branch_with_pr: true -# skip unsupported combinations init: - echo %APPVEYOR_BUILD_WORKER_IMAGE% -- if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2013" ( set generator="Visual Studio 12 2013" ) -- if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2015" ( set generator="Visual Studio 14 2015" ) - if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2017" ( set generator="Visual Studio 15 2017" ) +- if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2022" ( set generator="Visual Studio 17 2022" ) - echo %generator% before_build: diff --git a/src/include/cppduals/.gitlab-ci.yml b/src/include/cppduals/.gitlab-ci.yml index 9296366be..8dac950f8 100644 --- a/src/include/cppduals/.gitlab-ci.yml +++ b/src/include/cppduals/.gitlab-ci.yml @@ -1,5 +1,4 @@ -#image: ubuntu:19.04 -image: fedora:30 +image: fedora:41 variables: GIT_DEPTH: 3 @@ -9,79 +8,105 @@ stages: - build - test - cover -- publish - -before_script: -#- apt-get update --yes -#- apt-get install --yes cmake g++ git doxygen lcov graphviz -- dnf install -y gcc-c++ make cmake git doxygen lcov graphviz +- pages +- upload +- release build: stage: build -# variables: -# CC: clang -# CXX: clang++ script: + - dnf install -y gcc-c++ make cmake git doxygen lcov graphviz - echo $CXX - - cmake -Bbuild -H. -DCPPDUALS_TESTING=ON + - gcc -march=native -dM -E - // SSE3 namespace Eigen { @@ -150,17 +151,9 @@ template<> EIGEN_STRONG_INLINE void prefetch >(const duals::d template<> EIGEN_STRONG_INLINE duals::dual pfirst(const Packet2df& a) { - #if EIGEN_GNUC_AT_MOST(4,3) - // Workaround gcc 4.2 ICE - this is not performance wise ideal, but who cares... - // This workaround also fix invalid code generation with gcc 4.3 - EIGEN_ALIGN16 duals::dual res[2]; - _mm_store_ps((float*)res, a.v); - return res[0]; - #else duals::dual res; _mm_storel_pi((__m64*)&res, a.v); return res; - #endif } template<> EIGEN_STRONG_INLINE Packet2df preverse(const Packet2df& a) diff --git a/src/include/cppduals/duals/dual b/src/include/cppduals/duals/dual index 4f5bf8acc..94688d185 100644 --- a/src/include/cppduals/duals/dual +++ b/src/include/cppduals/duals/dual @@ -16,14 +16,24 @@ #define CPPDUALS_DUAL #ifndef PARSED_BY_DOXYGEN -#include -#include +#define _DUALS_CONSTEXPR constexpr + +///// FOR LIBCPP ///// +#include +#if !defined(_LIBCPP_BEGIN_NAMESPACE_STD) +namespace std { +template class complex; +} +#else // LIBCPP +_LIBCPP_BEGIN_NAMESPACE_STD +template class complex; +_LIBCPP_END_NAMESPACE_STD +#endif // LIBCPP + #include #include -#include -#include -#include -#endif + +#endif // PARSED_BY_DOXYGEN #if !defined(CPPDUALS_IGNORE_COMPILER_VERSION) && !defined(_WIN32) #if __cplusplus < 201103L @@ -34,11 +44,9 @@ /// 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 { /** @@ -180,11 +188,12 @@ T rpart(dual) // Real part T dpart(dual) // Dual part dual dconj(dual) // Dual-conjugate -dual random(dual low = {0,0}, dual high = {1,0}) - +// Transcendental functions dual exp(dual) dual log(dual) dual log10(dual) +dual log2(dual) +dual logb(dual x) dual pow(dual, U) dual pow(U, dual) dual pow(dual, dual) @@ -196,7 +205,9 @@ dual tan(dual) dual asin(dual) dual acos(dual) dual atan(dual) - +dual atan2(dual, dual) +dual hypot(dual, dual) +dual scalbn(dual arg, int ex) // TODO: dual sinh(dual) dual cosh(dual) @@ -205,18 +216,35 @@ dual asinh(dual) dual acosh(dual) dual atanh(dual) +dual erf(dual) +dual erfc(dual) +dual tgamma(dual) +dual lgamma(dual) + // Non-differentiable operations on the real part. -T frexp(duals::dual arg, int* exp ); -duals::dual ldexp(duals::dual arg, int exp ); -T trunc(duals::dual d); -T floor(duals::dual d); -T ceil(duals::dual d); -T round(duals::dual d); -int fpclassify(duals::dual d); -bool isfinite(duals::dual d); -bool isnormal(duals::dual d); -bool isinf(duals::dual d); -bool isnan(duals::dual d); +dual abs(dual d) +dual fabs(dual d) +dual fmax(dual x, dual y) +dual fmin(dual x, dual y) +dual frexp(dual arg, int* exp ) +dual ldexp(dual arg, int exp ) +dual trunc(dual d) +dual floor(dual d) +dual ceil(dual d) +dual round(dual d) + +// floating point functions +int fpclassify(dual d) +bool isfinite(dual d) +bool isnormal(dual d) +bool isinf(dual d) +bool isnan(dual d) +bool signbit(dual d) +dual copysign(dual x, dual y) + +// Utility functions +dual random(dual low = {0,0}, dual high = {1,0}) +dual random2(dual low = {0,0}, dual high = {1,0}) // Stream IO template operator>>(basic_istream &, dual &); @@ -256,7 +284,7 @@ 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: +printed, these must come before the normal `{fmt}` formatting spec '$', '*', ','. For me these are about 3x faster than iostreams. ``` @@ -333,7 +361,7 @@ template struct dual_traits> template struct dual_traits>> { /// complex> have the same 'depth' as their dual. - enum { depth = dual_traits::depth }; + enum { depth = dual_traits >::depth }; /// The real storage type. typedef typename dual_traits::real_type real_type; @@ -361,6 +389,8 @@ template struct promote : std::common_type using can_promote = detail::has_member_type>; +// dual dual +// T == U template struct promote, dual, typename std::enable_if<(can_promote::value @@ -368,6 +398,9 @@ struct promote, dual, { typedef dual::type> type; }; +// dual U +// T >= U +// U is not complex template struct promote, U, typename std::enable_if<(can_promote::value @@ -376,6 +409,9 @@ struct promote, U, { typedef dual::type> type; }; +// U dual +// T >= U +// U is not complex template struct promote, typename std::enable_if<(can_promote::value @@ -385,6 +421,8 @@ struct promote, typedef dual::type> type; }; // ///////////////////////////////////////////////// +// complex complex +// T == U template struct promote, std::complex, typename std::enable_if<(can_promote::value @@ -392,6 +430,9 @@ struct promote, std::complex, { typedef std::complex::type> type; }; +// complex U +// T >= U +// U is not complex template struct promote, U, typename std::enable_if<(can_promote::value @@ -400,6 +441,9 @@ struct promote, U, { typedef std::complex::type> type; }; +// U complex +// T >= U +// U is not complex template struct promote, typename std::enable_if<(can_promote::value @@ -446,7 +490,7 @@ struct is_compound> : true_type {}; // Modification of std::numeric_limits<> per // C++03 17.4.3.1/1, and C++11 18.3.2.3/1. template -struct numeric_limits> : numeric_limits { +struct numeric_limits > : public numeric_limits { static constexpr bool is_specialized = true; static constexpr duals::dual min NOMACRO () { return numeric_limits::min NOMACRO (); } static constexpr duals::dual lowest() { return numeric_limits::lowest(); } @@ -477,24 +521,24 @@ namespace duals { // Both T and U are wrapped in a dual<> #define CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U) \ - typename std::enable_if<(int)duals::dual_traits::depth == \ - (int)duals::dual_traits::depth, int>::type = 0, \ - typename std::enable_if::value,int>::type = 0, \ - typename common_t = dual::type> + typename common_t = dual::type>, \ + typename std::enable_if<(int)duals::dual_traits::depth == \ + (int)duals::dual_traits::depth, int>::type = 0, \ + typename std::enable_if::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, U>::type, \ typename std::enable_if<((int)duals::dual_traits::depth >= \ (int)duals::dual_traits::depth), int>::type = 0, \ - typename std::enable_if,U>::value,int>::type = 0, \ - typename common_t = typename duals::promote, U>::type + typename std::enable_if,U>::value,int>::type = 0 // T is wrapped in complex> #define CPPDUALS_ENABLE_LEQ_DEPTH_AND_CX_COMMON_T(T,U) \ + typename common_t = typename duals::promote >,U>::type, \ typename std::enable_if<((int)duals::dual_traits::depth >= \ (int)duals::dual_traits::depth), int>::type = 0, \ - typename std::enable_if>,U>::value,int>::type = 0, \ - typename common_t = typename duals::promote >,U>::type + typename std::enable_if>,U>::value,int>::type = 0 #define CPPDUALS_ENABLE_IF(...) typename std::enable_if< (__VA_ARGS__) , int>::type = 0 @@ -528,16 +572,13 @@ public: : _real(re), _dual(du) {} /// Copy construct from a dual of equal depth. -#if defined (_MSC_VER) -#pragma warning (disable : 4244) /* floats are not used in HICUM */ -#endif template::value)> dual(const dual & x) : _real(x.rpart()), _dual(x.dpart()) {} /// Cast to a complex> with real part equal to *this. - operator std::complex>() { return std::complex>(*this); } + operator std::complex>() { return std::complex>(*this, (T)0); } /// Explicit cast to an arithmetic type retains the rpart() template dual & operator-=(const dual & x) { _real -= x.rpart(); _dual -= x.dpart(); return *this; } - /// Multiply this dual with a dual of same of lower depth. + /// Multiply this dual with a dual of same depth. template dual & operator*=(const dual & x) { _dual = _real * x.dpart() + _dual * x.rpart(); @@ -733,6 +774,768 @@ CPPDUALS_LHS_COMPARISON(!=) #endif // PARSED_BY_DOXYGEN +/// Complex Conjugate of a dual is just the dual. +template dual conj(const dual & 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. +template::value && + !is_complex::value && + std::is_arithmetic::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 dual dconj(const dual & x) { + return dual(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 std::complex dconj(const std::complex & x) { + return std::complex(dconj(x.real()), dconj(x.imag())); +} + +/// Conjugate a thing that's not dual and not complex. +template::value && + !is_complex::value && + std::is_arithmetic::value)> +T dconj(const T & x) { return x; } + +namespace detail { + template dual exp_impl(const dual & x); + template dual log_impl(const dual & x); + template dual log10_impl(const dual & x); + template dual log2_impl(const dual & x); + template dual logb_impl(const dual & x); + template + common_t pow_dual_scalar(const dual& x, const U& y); + template + common_t pow_scalar_dual(const U& x, const dual& y); + template + std::complex> pow_complex_dual(const std::complex>& x, const dual& y); + template + std::complex> pow_dual_complex(const dual& x, const std::complex>& y); + template + common_t pow_dual_dual(const dual& x, const dual& y); + template dual sqrt_impl(const dual & x); + template dual cbrt_impl(const dual & x); + template dual sin_impl(const dual & x); + template dual cos_impl(const dual & x); + template dual tan_impl(const dual & x); + template dual asin_impl(const dual & x); + template dual acos_impl(const dual & x); + template dual atan_impl(const dual & x); + template + common_t atan2_dual_scalar(const dual & y, const U & x); + template + common_t atan2_scalar_dual(const U & y, const dual & x); + template dual atan2_dual_dual(const dual & y, const dual & x); + template dual hypot_impl(const dual & x, const dual & y); + template dual scalbn_impl(const dual & arg, int ex); + template dual sinh_impl(const dual & x); + template dual cosh_impl(const dual & x); + template dual tanh_impl(const dual & x); + template dual asinh_impl(const dual & x); + template dual acosh_impl(const dual & x); + template dual atanh_impl(const dual & x); + template dual erf_impl(const dual & x); + template dual erfc_impl(const dual & x); + template dual tgamma_impl(const dual & x); + template dual lgamma_impl(const dual & x); + template dual abs_impl(const dual & x); + template dual fabs_impl(const dual & x); + template dual fmax_impl(const dual & x, const dual & y); + template dual fmin_impl(const dual & x, const dual & y); + template dual frexp_impl(const dual & x, int* exp); + template dual ldexp_impl(const dual & x, int exp); + template dual trunc_impl(const dual & x); + template dual floor_impl(const dual & x); + template dual ceil_impl(const dual & x); + template dual round_impl(const dual & x); + template int fpclassify_impl(const dual & x); + template bool isfinite_impl(const dual & x); + template bool isnormal_impl(const dual & x); + template bool isinf_impl(const dual & x); + template bool isnan_impl(const dual & x); + template bool signbit_impl(const dual & x); + template dual copysign_impl(const dual & x, const dual & y); + template dual log1p_impl(const dual & x); + template dual expm1_impl(const dual & x); + template dual random_impl(const dual & low, const dual & high); + template dual random2_impl(const dual & low, const dual & high); +} + +/// /// /// /// /// /// /// /// +/// Prototypes +/// + +// Public interfaces that forward to implementations +template dual exp(const dual & x) { return detail::exp_impl(x); } +template dual log(const dual & x) { return detail::log_impl(x); } +template dual log10(const dual & x) { return detail::log10_impl(x); } +template dual log2(const dual & x) { return detail::log2_impl(x); } +template dual logb(const dual & x) { return detail::logb_impl(x); } +template::value)> +common_t pow(const dual& x, const U& y) { return detail::pow_dual_scalar(x, y); } +template::value)> +common_t pow(const U& x, const dual& y) { return detail::pow_scalar_dual(x, y); } +template +std::complex> pow(const dual & x, const std::complex> & y) { return detail::pow_dual_complex(x, y); } +template +std::complex> pow(const std::complex> & x, const dual& y) { return detail::pow_complex_dual(x, y); } +template +common_t pow(const dual& x, const dual& y) { return detail::pow_dual_dual(x, y); } +template dual sqrt(const dual & x) { return detail::sqrt_impl(x); } +template dual cbrt(const dual & x) { return detail::cbrt_impl(x); } +template dual sin(const dual & x) { return detail::sin_impl(x); } +template dual cos(const dual & x) { return detail::cos_impl(x); } +template dual tan(const dual & x) { return detail::tan_impl(x); } +template dual asin(const dual & x) { return detail::asin_impl(x); } +template dual acos(const dual & x) { return detail::acos_impl(x); } +template dual atan(const dual & x) { return detail::atan_impl(x); } +template +common_t atan2(const dual & y, const U & x) { return detail::atan2_dual_scalar(y, x); } +template +common_t atan2(const U & y, const dual & x) { return detail::atan2_scalar_dual(y, x); } +template dual atan2(const dual & y, const dual & x) { return detail::atan2_dual_dual(y, x); } +template dual hypot(const dual & x, const dual & y) { return detail::hypot_impl(x, y); } +template dual scalbn(const dual & arg, int ex) { return detail::scalbn_impl(arg, ex); } +template dual sinh(const dual & x) { return detail::sinh_impl(x); } +template dual cosh(const dual & x) { return detail::cosh_impl(x); } +template dual tanh(const dual & x) { return detail::tanh_impl(x); } +template dual asinh(const dual & x) { return detail::asinh_impl(x); } +template dual acosh(const dual & x) { return detail::acosh_impl(x); } +template dual atanh(const dual & x) { return detail::atanh_impl(x); } +/// The error function. Make sure to `#include ` before +/// `#include ` to use this function. +template dual erf(const dual & x) { return detail::erf_impl(x); } +template dual erfc(const dual & x) { return detail::erfc_impl(x); } +template dual tgamma(const dual & x) { return detail::tgamma_impl(x); } +template dual lgamma(const dual & x) { return detail::lgamma_impl(x); } +template dual abs(const dual & x) { return detail::abs_impl(x); } +template dual fabs(const dual & x) { return detail::fabs_impl(x); } +template dual fmax(const dual & x, const dual & y) { return detail::fmax_impl(x, y); } +template dual fmin(const dual & x, const dual & y) { return detail::fmin_impl(x, y); } + +template dual frexp(const dual & x, int* exp) { return detail::frexp_impl(x, exp); } +template dual ldexp(const dual & x, int exp) { return detail::ldexp_impl(x, exp); } +//template dual trunc(const dual & x) { return detail::trunc_impl(x); } +template dual trunc(const dual & x) { return detail::trunc_impl(x); } +template dual floor(const dual & x) { return detail::floor_impl(x); } +template dual ceil(const dual & x) { return detail::ceil_impl(x); } +template dual round(const dual & x) { return detail::round_impl(x); } + +template int fpclassify(const dual & x) { return detail::fpclassify_impl(x); } +template bool isfinite(const dual & x) { return detail::isfinite_impl(x); } +template bool isnormal(const dual & x) { return detail::isnormal_impl(x); } +template bool isinf(const dual & x) { return detail::isinf_impl(x); } +template bool isnan(const dual & x) { return detail::isnan_impl(x); } +template bool signbit(const dual & x) { return detail::signbit_impl(x); } +template dual copysign(const dual & x, const dual & y) { return detail::copysign_impl(x, y); } + +//template dual random(const dual & low, const dual & high) { return detail::random_impl(low, high); } +//template dual random2(const dual & low, const dual & high) { return detail::random2_impl(low, high); } + +namespace utils { +template int sgn(T val) { + return (T(0) < val) - (val < T(0)); +} +} + +template dual log1p(const dual & x) {return detail::log1p_impl(x);} +template dual expm1(const dual & 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 operator""_ef(long double du) +{ + return { 0.0f, static_cast(du) }; +} +constexpr dual operator""_ef(unsigned long long du) +{ + return { 0.0f, static_cast(du) }; +} +constexpr dual operator""_e(long double du) +{ + return { 0.0, static_cast(du) }; +} +constexpr dual operator""_e(unsigned long long du) +{ + return { 0.0, static_cast(du) }; +} +constexpr dual operator""_el(long double du) +{ + return { 0.0l, du }; +} +constexpr dual operator""_el(unsigned long long du) +{ + return { 0.0l, static_cast(du) }; +} +} // namespace literals +#endif + +// shortcut type names +typedef dual dualf; +typedef dual duald; +typedef dual dualld; +typedef dual hyperdualf; +typedef dual hyperduald; +typedef dual 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 _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_copysign(duals::dual
a, duals::dual
b); +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_fabs(duals::dual
a); +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isfinite(duals::dual
a); +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isinf(duals::dual
a); +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_fmax(duals::dual
a, duals::dual
b); +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_scalbn(duals::dual
a, int exp); +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_logb(duals::dual
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 +#include + +#ifdef _LIBCPP_BEGIN_NAMESPACE_STD +_LIBCPP_BEGIN_NAMESPACE_STD +#else +namespace std { +#endif + +#if defined(_LIBCPP_BEGIN_NAMESPACE_STD) +// special needs LIBCPP +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_copysign(duals::dual
a, duals::dual
b) { return duals::copysign(a, b); } +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_fabs(duals::dual
a) { return duals::fabs(a); } +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isfinite(duals::dual
a) { return duals::isfinite(a); } +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI bool __constexpr_isinf(duals::dual
a) { return duals::isinf(a); } +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_fmax(duals::dual
a, duals::dual
b) { return duals::fmax(a, b); } +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_scalbn(duals::dual
a, int exp) { return duals::scalbn(a, exp); } +template _LIBCPP_CONSTEXPR inline _LIBCPP_HIDE_FROM_ABI duals::dual
__constexpr_logb(duals::dual
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 dual exp_impl(const dual & x) { + using std::exp; + T v = exp(x.rpart()); + return dual(v, + v * x.dpart()); +} + +/// Natural log ln(x) +template dual log_impl(const dual & x) { + using std::log; + T v = log(x.rpart()); + if (x.dpart() == T(0)) + return v; + else + return dual(v, x.dpart() / x.rpart()); +} + +template dual log10_impl(const dual & x) { + using std::log; + return log(x) / log(static_cast(10)); +} + +template dual log2_impl(const dual & x) { + using std::log; + return log(x) / log(static_cast(2)); +} + +template duals::dual logb_impl(const duals::dual & x) { + using std::log; + using std::pow; + using std::abs; + using std::floor; + T v = floor(log(abs(x.rpart())) / log(std::numeric_limits::radix)); + return duals::dual(v, + x.rpart() - pow(2,v) ? (T)0 : std::numeric_limits::infinity()); +} + +// dual ^ U implementation +template +common_t pow_dual_scalar(const dual& 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 implementation +template +common_t pow_scalar_dual(const U& x, const dual& y) +{ + return pow(common_t(x), y); +} + +// U ^ dual implementation +template +std::complex> pow_complex_dual(const std::complex>& x, const dual& y) +{ + using dual_t = duals::dual; + + // 1. Compute the magnitude (r = |x|). + // r = sqrt( (Re x)^2 + (Im x)^2 ), where Re x and Im x are dual. + 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>. + return std::complex(realOut, imagOut); +} + +template +std::complex> pow_dual_complex( + const dual& x, const std::complex>& y) +{ + // For clarity, define short aliases: + using dual_t = duals::dual; + + // Extract the real and imaginary parts of y: both are dual. + 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(realPart, imagPart); +} + +// dual ^ dual implementation +template +common_t pow_dual_dual(const dual& f, const dual& 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::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 dual sqrt_impl(const dual & x) { + using std::sqrt; + T v = sqrt(x.rpart()); + if (x.dpart() == T(0)) + return v; + else + return dual(v, x.dpart() / (T(2) * v) ); +} + +template dual cbrt_impl(const dual & x) { + using std::cbrt; + T v = cbrt(x.rpart()); + if (x.dpart() == T(0)) + return v; + else + return dual(v, x.dpart() / (T(3) * v * v) ); +} + +template dual sin_impl(const dual & x) { + using std::sin; + using std::cos; + return dual(sin(x.rpart()), + x.dpart() * cos(x.rpart())); +} + +template dual cos_impl(const dual & x) { + using std::cos; + using std::sin; + return dual(cos(x.rpart()), + -sin(x.rpart()) * x.dpart()); +} + +template dual tan_impl(const dual & x) { + using std::tan; + T v = tan(x.rpart()); + return dual(v, x.dpart() * (v*v + 1)); +} + +template dual asin_impl(const dual & x) { + using std::asin; + using std::sqrt; + T v = asin(x.rpart()); + if (x.dpart() == T(0)) + return v; + else + return dual(v, x.dpart() / sqrt(1 - x.rpart()*x.rpart())); +} + +template dual acos_impl(const dual & x) { + using std::acos; + using std::sqrt; + T v = acos(x.rpart()); + if (x.dpart() == T(0)) + return v; + else + return dual(v, -x.dpart() / sqrt(1 - x.rpart()*x.rpart())); +} + +template dual atan_impl(const dual & x) { + using std::atan; + T v = atan(x.rpart()); + if (x.dpart() == T(0)) + return v; + else + return dual(v, x.dpart() / (1 + x.rpart() * x.rpart())); +} + +template dual atan2_dual_dual(const dual & y, const dual & x) { + using std::atan2; + T v = atan2(y.rpart(), x.rpart()); + return dual(v, + (x.rpart() * y.dpart() - y.rpart() * x.dpart()) / + (y.rpart() * y.rpart() + x.rpart() * x.rpart())); +} + +template +common_t +atan2_dual_scalar(const dual & y, const U & x) { + using std::atan2; + T v = atan2(y.rpart(), x); + return dual(v, + (x * y.dpart()) / + (y.rpart() * y.rpart() + x * x)); +} + +template +common_t +atan2_scalar_dual(const U & y, const dual & x) { + using std::atan2; + T v = atan2(y, x.rpart()); + return dual(v, + (-y * x.dpart()) / + (y * y + x.rpart() * x.rpart())); +} + +template duals::dual hypot_impl(const duals::dual & x, const duals::dual & y) { + using std::hypot; + return dual(hypot(x.rpart(), y.rpart()), + (x.rpart() * y.dpart() + y.rpart() * x.dpart()) / hypot(x.rpart(), y.rpart())); +} + +template duals::dual scalbn_impl(const duals::dual & arg, int ex) { + return arg * std::pow((T)2, ex); +} + +template dual sinh_impl(const dual & x) { + using std::sinh; + using std::cosh; + return dual(sinh(x.rpart()), + x.dpart() * cosh(x.rpart())); +} + +template dual cosh_impl(const dual & x) { + using std::cosh; + using std::sinh; + return dual(cosh(x.rpart()), + x.dpart() * sinh(x.rpart())); +} + +template dual tanh_impl(const dual & x) { + using std::tanh; + return dual(tanh(x.rpart()), + x.dpart() * (1 - tanh(x.rpart()) * tanh(x.rpart()))); +} + +template dual asinh_impl(const dual & x) { + using std::asinh; + using std::sqrt; + return dual(asinh(x.rpart()), + x.dpart() / sqrt(1 + x.rpart() * x.rpart())); +} + +template dual acosh_impl(const dual & x) { + using std::acosh; + using std::sqrt; + return dual(acosh(x.rpart()), + x.dpart() / sqrt(x.rpart() * x.rpart() - 1)); +} + +template dual atanh_impl(const dual & x) { + using std::atanh; + return dual(atanh(x.rpart()), + x.dpart() / (1 - x.rpart() * x.rpart())); +} + +template dual log1p_impl(const dual & x) { + using std::log1p; + return dual(log1p(x.rpart()), + x.dpart() / (1 + x.rpart())); +} + +template dual expm1_impl(const dual & x) { + using std::expm1; + using std::exp; + return dual(expm1(x.rpart()), + x.dpart() * exp(x.rpart())); +} + +/// The error function. Make sure to `#include ` before +/// `#include ` to use this function. +template dual erf_impl(const dual & x) { + using std::erf; + using std::sqrt; + using std::pow; + using std::exp; + return dual(erf(x.rpart()), + x.dpart() * T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2)))); +} + +/// Error function complement (1 - erf()). +template dual erfc_impl(const dual & x) { + using std::erfc; + using std::sqrt; + using std::pow; + using std::exp; + return dual(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 dual tgamma_impl(const dual & 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::digits / 3))); + T w((tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h)); + errno = errno_saved; + return dual(v, x.dpart() * w); + } +} + +/// Log of absolute value of gamma function. Approximation of the dual part. +template dual lgamma_impl(const dual & 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::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(v, x.dpart() * w); + } +} +template dual abs_impl(const dual & x) { + using std::abs; + return dual(abs(x.rpart()), + x.dpart() * utils::sgn(x.rpart())); +} +template dual fabs_impl(const dual & x) { + using std::fabs; + return dual(fabs(x.rpart()), + x.dpart() * utils::sgn(x.rpart())); +} +template duals::dual fmax_impl(const duals::dual & x, const duals::dual & y) { + return x.rpart() > y.rpart() ? x : y; +} +template duals::dual fmin_impl(const duals::dual & x, const duals::dual & y) { + return x.rpart() <= y.rpart() ? x : y; +} +template duals::dual frexp_impl(const duals::dual & arg, int* exp ) { using std::frexp; return frexp(arg.rpart(), exp); } +template duals::dual ldexp_impl(const duals::dual & arg, int exp ) { return arg * std::pow((T)2,exp); } +template duals::dual trunc_impl(const duals::dual & d) { using std::trunc; return trunc(d.rpart()); } +template dual floor_impl(const dual & x) { + using std::floor; + return dual(floor(x.rpart()), + x.dpart() * utils::sgn(x.rpart()-floor(x.rpart()))); +} +template duals::dual ceil_impl(const duals::dual & d) { using std::ceil; return ceil(d.rpart()); } +template duals::dual round_impl(const duals::dual & d) { using std::round; return round(d.rpart()); } +template int fpclassify_impl(const duals::dual & d) { using std::fpclassify; return fpclassify(d.rpart()); } +template bool isfinite_impl(const duals::dual & d) { using std::isfinite; return isfinite(d.rpart()); } +template bool isnormal_impl(const duals::dual & d) { using std::isnormal; return isnormal(d.rpart()); } +template bool isinf_impl(const duals::dual & d) { using std::isinf; return isinf(d.rpart()); } +template bool isnan_impl(const duals::dual & d) { using std::isnan; return isnan(d.rpart()); } +template bool signbit_impl(const duals::dual & d) { using std::signbit; return signbit(d.rpart()); } +template duals::dual copysign_impl(const duals::dual & x, const duals::dual & y) { + using std::copysign; + T r = copysign(x.rpart(), y.rpart()); + return duals::dual(r, r == x.rpart() ? x.dpart() : -x.dpart()); +} + +} // namespace detail +} // namespace duals + +#include +namespace duals { namespace randos { // Random real value between a and b. @@ -781,290 +1584,18 @@ DT random(DT a = DT(0,0), DT b = DT(1,0)) { random(a.dpart(),b.dpart())); } -/// Complex Conjugate of a dual is just the dual. -template dual conj(const dual & x) { return x; } +} // namespace duals + +#include -/// 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. -template::value && - !is_complex::value && - std::is_arithmetic::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 dual dconj(const dual & x) { - return dual(x.rpart(), - x.dpart()); +namespace duals { +typedef std::complex cdualf; +typedef std::complex cduald; +typedef std::complex cdualld; } -/// Dual Conjugate of a complex, such that x*dconj(x) = rpart(x)^2. -/// Different function name from complex conjugate conj(). -template std::complex dconj(const std::complex & x) { - return std::complex(dconj(x.real()), dconj(x.imag())); -} - -/// Conjugate a thing that's not dual and not complex. -template::value && - !is_complex::value && - std::is_arithmetic::value)> -T dconj(const T & x) { return x; } - -/// Exponential e^x -template dual exp(const dual & x) { - using std::exp; - T v = exp(x.rpart()); - return dual(v, - v * x.dpart()); -} - -/// Natural log ln(x) -template dual log(const dual & x) { - using std::log; - T v = log(x.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, x.dpart() / x.rpart()); -} - -template dual log10(const dual & x) { - using std::log; - return log(x) / log(static_cast(10)); -} - -template dual log2(const dual & x) { - using std::log; - return log(x) / log(static_cast(2)); -} - -template -common_t -pow(const dual & x, const dual & y) { - using std::pow; - using std::log; - T v = pow(x.rpart(), y.rpart()); - return common_t(v, - x.dpart() * y.rpart() * pow(x.rpart(), y.rpart() - T(1)) + - y.dpart() * log(x.rpart()) * v); -} - -template -common_t -pow(const dual & x, const U & y) { - using std::pow; - return common_t(pow(x.rpart(), y), - x.dpart() * y * pow(x.rpart(), y - U(1))); -} - -template -common_t -pow(const U & x, const dual & y) { - return pow(common_t(x), y); -} - -namespace utils { -template int sgn(T val) { - return (T(0) < val) - (val < T(0)); -} -} - -template dual abs(const dual & x) { - using std::abs; - return dual(abs(x.rpart()), - x.dpart() * utils::sgn(x.rpart())); -} - -template dual fabs(const dual & x) { - using std::fabs; - return dual(fabs(x.rpart()), - x.dpart() * utils::sgn(x.rpart())); -} - -#if 0 -template dual abs2(const dual & x) { - using std::abs; - return dual(x.rpart() * x.rpart(), - xxx x.dpart() * utils::sgn(x.rpart())); -} -#endif - -template duals::dual copysign(const duals::dual & x, const duals::dual & y) { - using std::copysign; - T r = copysign(x.rpart(), y.rpart()); - return duals::dual(r, r == x.rpart() ? x.dpart() : -x.dpart()); -} - -template duals::dual hypot(const duals::dual & x, const duals::dual & y) { - return sqrt(x*x + y*y); -} - -template duals::dual scalbn(const duals::dual & arg, int ex) { - return arg * std::pow((T)2, ex); -} - -template duals::dual (fmax)(const duals::dual & x, const duals::dual & y) { - return x.rpart() > y.rpart() ? x : y; -} - -template duals::dual (fmin)(const duals::dual & x, const duals::dual & y) { - return x.rpart() <= y.rpart() ? x : y; -} - -template duals::dual logb(const duals::dual & x) { - return duals::log2(x); -} - -template int (fpclassify)(const duals::dual & d) { using std::fpclassify; return (fpclassify)(d.rpart()); } -template bool (isfinite)(const duals::dual & d) { using std::isfinite; return (isfinite)(d.rpart()); } -template bool (isnormal)(const duals::dual & d) { using std::isnormal; return (isnormal)(d.rpart()); } -template bool (isinf)(const duals::dual & d) { using std::isinf; return (isinf)(d.rpart()); } -template bool (isnan)(const duals::dual & d) { using std::isnan; return (isnan)(d.rpart()); } -template bool (signbit)(const duals::dual & d) { using std::signbit; return (signbit)(d.rpart()); } - -template dual sqrt(const dual & x) { - using std::sqrt; - T v = sqrt(x.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, x.dpart() / (T(2) * v) ); -} - -template dual cbrt(const dual & x) { - using std::cbrt; - T v = cbrt(x.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, x.dpart() / (T(3) * v * v) ); -} - -template dual sin(const dual & x) { - using std::sin; - using std::cos; - return dual(sin(x.rpart()), - x.dpart() * cos(x.rpart())); -} - -template dual cos(const dual & x) { - using std::cos; - using std::sin; - return dual(cos(x.rpart()), - -sin(x.rpart()) * x.dpart()); -} - -template dual tan(const dual & x) { - using std::tan; - T v = tan(x.rpart()); - return dual(v, x.dpart() * (v*v + 1)); -} - -template dual asin(const dual & x) { - using std::asin; - using std::sqrt; - T v = asin(x.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, x.dpart() / sqrt(1 - x.rpart()*x.rpart())); -} - -template dual acos(const dual & x) { - using std::acos; - using std::sqrt; - T v = acos(x.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, -x.dpart() / sqrt(1 - x.rpart()*x.rpart())); -} - -template dual atan(const dual & x) { - using std::atan; - T v = atan(x.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, x.dpart() / (1 + x.rpart()*x.rpart())); -} - -template dual atan2(const dual & x, const dual & y) { - using std::atan2; - T v = atan2(x.rpart(), y.rpart()); - if (x.dpart() == T(0)) - return v; - else - return dual(v, x.dpart() / (1 + x.rpart()*x.rpart())); -} - -// TODO -template dual sinh(const dual & x); -template dual cosh(const dual & x); -template dual tanh(const dual & x); -template dual asinh(const dual & x); -template dual acosh(const dual & x); -template dual atanh(const dual & x); -template dual log1p(const dual & x); -template dual expm1(const dual & x); - -/// The error function. Make sure to `#include ` before -/// `#include ` to use this function. -template dual erf(const dual & x) { - using std::erf; - using std::sqrt; - using std::pow; - using std::exp; - return dual(erf(x.rpart()), - x.dpart() * T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2)))); -} - -/// Error function complement (1 - erf()). -template dual erfc(const dual & x) { - using std::erfc; - using std::sqrt; - using std::pow; - using std::exp; - return dual(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 dual tgamma(const dual & 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::digits / 3))); - T w((tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h)); - errno = errno_saved; - return dual(v, x.dpart() * w); - } -} - -/// Log of absolute value of gamma function. Approximation of the dual part. -template dual lgamma(const dual & 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::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(v, x.dpart() * w); - } -} +#include +namespace duals { /// Putto operator template @@ -1170,48 +1701,6 @@ operator>>(std::basic_istream & is, dual & x) return is; } -#if __cpp_user_defined_literals >= 200809 -/// Dual number literals in namespace duals::literals -inline namespace literals -{ -using duals::dual; -constexpr dual operator""_ef(long double du) -{ - return { 0.0f, static_cast(du) }; -} -constexpr dual operator""_ef(unsigned long long du) -{ - return { 0.0f, static_cast(du) }; -} -constexpr dual operator""_e(long double du) -{ - return { 0.0, static_cast(du) }; -} -constexpr dual operator""_e(unsigned long long du) -{ - return { 0.0, static_cast(du) }; -} -constexpr dual operator""_el(long double du) -{ - return { 0.0l, du }; -} -constexpr dual operator""_el(unsigned long long du) -{ - return { 0.0l, static_cast(du) }; -} -} -#endif - -typedef dual dualf; -typedef dual duald; -typedef dual dualld; -typedef dual hyperdualf; -typedef dual hyperduald; -typedef dual hyperdualld; -typedef std::complex cdualf; -typedef std::complex cduald; -typedef std::complex cdualld; - } // namespace duals #ifdef CPPDUALS_LIBFMT @@ -1221,7 +1710,7 @@ typedef std::complex cdualld; /// 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 +/// 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 @@ -1236,47 +1725,61 @@ struct fmt::formatter,Char> : public fmt::formatter { typedef fmt::formatter base; enum style { expr, star, pair } style_ = expr; - internal::dynamic_format_specs specs_; - FMT_CONSTEXPR auto parse(format_parse_context & ctx) -> decltype(ctx.begin()) { - using handler_type = internal::dynamic_specs_handler; - auto type = internal::type_constant::value; - internal::specs_checker handler(handler_type(specs_, ctx), type); + detail::dynamic_format_specs specs_; + + FMT_CONSTEXPR auto parse(parse_context& ctx) -> decltype(ctx.begin()) { auto it = ctx.begin(); - 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; + 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; + } } - parse_format_specs(ctx.begin(), ctx.end(), handler); - return base::parse(ctx); + return parse_format_specs(it, ctx.end(), specs_, ctx, + detail::type_constant::value); } + template - auto format(const duals::dual & x, FormatCtx & ctx) -> decltype(ctx.out()) { - format_to(ctx.out(), "("); + auto format(const duals::dual & 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) { - base::format(x.rpart(), ctx); - format_to(ctx.out(), ","); - base::format(x.dpart(), ctx); - return format_to(ctx.out(), ")"); + out = detail::write(out, x.rpart(), specs, ctx.locale()); + *out++ = ','; + out = detail::write(out, x.dpart(), specs, ctx.locale()); + } else { + if (x.rpart() || !x.dpart()) + out = detail::write(out, x.rpart(), specs, ctx.locale()); + if (x.dpart()) { + if (x.rpart() && x.dpart() >= 0) + specs.set_sign(sign::plus); + out = detail::write(out, x.dpart(), specs, ctx.locale()); + if (style_ == style::star) { + *out++ = '*'; + *out++ = 'e'; + } else { + *out++ = '_'; + *out++ = 'e'; + } + if (std::is_same::type,float>::value) *out++ = 'f'; + if (std::is_same::type,long double>::value) *out++ = 'l'; + } } - if (x.rpart() || !x.dpart()) - base::format(x.rpart(), ctx); - if (x.dpart()) { - if (x.rpart() && x.dpart() >= 0 && specs_.sign != sign::plus) - format_to(ctx.out(), "+"); - base::format(x.dpart(), ctx); - if (style_ == style::star) - format_to(ctx.out(), "*e"); - else - format_to(ctx.out(), "_e"); - if (std::is_same::type,float>::value) format_to(ctx.out(), "f"); - if (std::is_same::type,long double>::value) format_to(ctx.out(), "l"); - } - return format_to(ctx.out(), ")"); + *out++ = ')'; + return out; } }; -#endif +#endif // CPPDUALS_LIBFMT #ifdef CPPDUALS_LIBFMT_COMPLEX #ifndef CPPDUALS_LIBFMT @@ -1294,7 +1797,7 @@ struct fmt::formatter,Char> : public fmt::formatter /// 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 +/// (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). @@ -1308,82 +1811,53 @@ struct fmt::formatter,Char> : public fmt::formatter { typedef fmt::formatter base; enum style { expr, star, pair } style_ = expr; - internal::dynamic_format_specs specs_; - FMT_CONSTEXPR auto parse(format_parse_context & ctx) -> decltype(ctx.begin()) { - using handler_type = internal::dynamic_specs_handler; - auto type = internal::type_constant::value; - internal::specs_checker handler(handler_type(specs_, ctx), type); + detail::dynamic_format_specs specs_; + + FMT_CONSTEXPR auto parse(parse_context& ctx) -> decltype(ctx.begin()) { auto it = ctx.begin(); - 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; + 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; + } } - parse_format_specs(ctx.begin(), ctx.end(), handler); - //todo: fixup alignment - return base::parse(ctx); + // 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::value); + return end; } + template - auto format(const std::complex & x, FormatCtx & ctx) -> decltype(ctx.out()) { + auto format(const std::complex & 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); - return format_to(ctx.out(), ")"); - } - if (x.real() || !x.imag()) - base::format(x.real(), ctx); - if (x.imag()) { - if (x.real() && x.imag() >= 0 && specs_.sign != sign::plus) - 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::type,float>::value) format_to(ctx.out(), "f"); - if (std::is_same::type,long double>::value) format_to(ctx.out(), "l"); + } 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::type,float>::value) format_to(ctx.out(), "f"); + if (std::is_same::type,long double>::value) format_to(ctx.out(), "l"); + } } return format_to(ctx.out(), ")"); } }; -#endif - -#ifdef _LIBCPP_BEGIN_NAMESPACE_STD -_LIBCPP_BEGIN_NAMESPACE_STD -#else -namespace std { -#endif - -#ifndef PARSED_BY_DOXYGEN - -#define make_math(T) \ - inline T (frexp)(const duals::dual & arg, int* exp ) { return (frexp)(arg.rpart(), exp); } \ - inline duals::dual (ldexp)(const duals::dual & arg, int exp ) { return arg * std::pow((T)2,exp); } \ - inline T (trunc)(const duals::dual & d) { return (trunc)(d.rpart()); } \ - inline T (floor)(const duals::dual & d) { return (floor)(d.rpart()); } \ - inline T (ceil)(const duals::dual & d) { return (ceil)(d.rpart()); } \ - inline T (round)(const duals::dual & d) { return (round)(d.rpart()); } \ - inline int (fpclassify)(const duals::dual & d) { return (fpclassify)(d.rpart()); } \ - inline bool (isfinite)(const duals::dual & d) { return (isfinite)(d.rpart()); } \ - inline bool (isnormal)(const duals::dual & d) { return (isnormal)(d.rpart()); } \ - inline bool (isinf)(const duals::dual & d) { return (isinf)(d.rpart()); } \ - inline bool (isnan)(const duals::dual & d) { return (isnan)(d.rpart()); } \ - -make_math(float) -make_math(double) -make_math(long double) - -#undef make_math - -#endif // PARSED_BY_DOXYGEN - -#ifdef _LIBCPP_BEGIN_NAMESPACE_STD -_LIBCPP_END_NAMESPACE_STD -#else -} // namespace std -#endif +#endif // CPPDUALS_LIBFMT_COMPLEX #endif // CPPDUALS_DUAL diff --git a/src/include/cppduals/duals/dual_eigen b/src/include/cppduals/duals/dual_eigen index 818c8e7f7..cab406aa4 100644 --- a/src/include/cppduals/duals/dual_eigen +++ b/src/include/cppduals/duals/dual_eigen @@ -141,7 +141,6 @@ namespace Eigen { template struct NumTraits > : GenericNumTraits { - typedef typename NumTraits::Real ReallyReal; typedef duals::dual Real; typedef duals::dual Literal; typedef duals::dual Nested; @@ -159,13 +158,13 @@ struct NumTraits > : GenericNumTraits EIGEN_DEVICE_FUNC static inline Real epsilon() { return Real(NumTraits::epsilon()); } EIGEN_DEVICE_FUNC - static inline ReallyReal dummy_precision() { return NumTraits::dummy_precision(); } + static inline Real dummy_precision() { return NumTraits::dummy_precision(); } EIGEN_DEVICE_FUNC - static inline ReallyReal highest() { return NumTraits::highest(); } + static inline Real highest() { return NumTraits::highest(); } EIGEN_DEVICE_FUNC - static inline ReallyReal lowest() { return NumTraits::lowest(); } + static inline Real lowest() { return NumTraits::lowest(); } EIGEN_DEVICE_FUNC - static inline ReallyReal digits10() { return NumTraits::digits10(); } + static inline int digits10() { return NumTraits::digits10(); } }; #if !defined(CPPDUALS_NO_EIGEN_PROMOTION) @@ -218,6 +217,20 @@ using duals::dconj; namespace internal { +#if 0 +// For MatrixExponential.h to treat duals::dual as a known type and compile it +template struct is_exp_known_type; +template struct matrix_exp_computeUV; +template struct is_exp_known_type> : is_exp_known_type {}; +template +struct matrix_exp_computeUV > : matrix_exp_computeUV +{ + typedef typename NumTraits::Scalar>::Real RealScalar; +}; +#endif + +#if 1 +// this is used by the packet math for SSE to copy raw duals around. template struct real_impl > { @@ -250,7 +263,44 @@ struct real_ref_retval> { typedef typename NumTraits::Real & type; }; + +#else //// +template +struct real_impl > +{ + typedef duals::dual RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const duals::dual& x) + { + return x; + } +}; + +template +struct real_ref_impl> +{ + typedef typename NumTraits >::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar & run(duals::dual & x) + { + return reinterpret_cast(&x)[0]; + } + EIGEN_DEVICE_FUNC + static inline const RealScalar & run(const duals::dual & x) + { + return reinterpret_cast(&x)[0]; + } +}; + +template +struct real_ref_retval> +{ + typedef typename NumTraits>::Real & type; +}; + +#endif /// 0 + template struct scalar_random_op> { EIGEN_EMPTY_STRUCT_CTOR(scalar_random_op) @@ -384,7 +434,8 @@ template struct dconj_helper\\n\", \"\", \"\", \"[\", \"]\")") -add_definitions (-DEIGEN_DEFAULT_IO_FORMAT=${IOFORMAT}) -#add_definitions (-DEIGEN_DEFAULT_IO_FORMAT=EIGEN_IO_FORMAT_OCTAVE) # # Correctness & Coverage # -if (NOT MSVC) - set (OPT_FLAGS "-O2") - set (OPT_FLAGS "-O3;-msse3") - set (OPT_FLAGS "-O3;-mavx2;-mfma") - set (OPT_FLAGS "-O3;-march=native") +message ("OSX_ARCHITECTURES: ${CMAKE_SYSTEM_PROCESSOR}") +set (OSX_ARCHITECTURES "arm64;x86_64") +if (CMAKE_SYSTEM_PROCESSOR MATCHES "(x86)|(X86)|(amd64)|(AMD64)") + set (X86 TRUE) else () - set (OPT_FLAGS "/arch:IA32") - set (OPT_FLAGS "/arch:SSE") - set (OPT_FLAGS "/arch:SSE2") - set (OPT_FLAGS "/arch:AVX2") + set (X86 FALSE) +endif () + +# Set optimization flags only if coverage is not enabled +if (NOT CPPDUALS_CODE_COVERAGE) + if (X86) + if (NOT MSVC) + set (OPT_FLAGS "-O3;-march=native") + else () + set (OPT_FLAGS "/arch:AVX2") + endif () + else () + set (OPT_FLAGS "-O3") + endif () +else() + if (X86) + if (NOT MSVC) + set (OPT_FLAGS "-mavx2") + else () + set (OPT_FLAGS "/arch:AVX2") + endif () + else () + set (OPT_FLAGS "") + endif () endif () -#set (OPT_FLAGS "${OPT_FLAGS};-fsanitize=address;-fno-omit-frame-pointer") set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_VECTORIZE_CDUAL") -#set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_DONT_VECTORIZE") -#set (OPT_FLAGS "${OPT_FLAGS};-DEIGEN_DONT_VECTORIZE") set (ALL_TESTS - test_dual test_funcs test_eigen test_packets + test_dual test_cdual test_funcs test_eigen test_packets test_vectorize test_solve test_expm test_1 test_fmt example ) @@ -70,33 +76,30 @@ foreach (TEST_ ${ALL_TESTS}) add_executable (${TEST} ${TEST_}.cpp) set (PHASE_FLAGS ${OPT_FLAGS} -DPHASE_${PHASE}) - string (REPLACE ";" ", " L2 "${CMAKE_CXX_FLAGS}") + string (REPLACE ";" ", " L2 "${OPT_FLAGS};${CMAKE_CXX_FLAGS}") target_compile_options (${TEST} PUBLIC ${PHASE_FLAGS}) target_compile_definitions (${TEST} PRIVATE "OPT_FLAGS=${L2}") - target_link_libraries (${TEST} gtest_main cppduals_coverage_config) - #target_link_libraries (${TEST} -lasan) - add_dependencies (${TEST} eigenX expokitX) + + # Link with coverage config first to ensure flags are used + target_link_libraries (${TEST} + cppduals_coverage_config # Link coverage first + gtest_main + eigen + expokit + ) + gtest_discover_tests (${TEST} TEST_LIST ${TEST}_targets) set_tests_properties (${${TEST}_targets} PROPERTIES TIMEOUT 10) - # -ftest-coverage endforeach (PHASE) endforeach (TEST_) # special for fmt -target_compile_features (test_fmt_1 PUBLIC cxx_std_14) target_link_libraries (test_fmt_1 fmt::fmt) if (CPPDUALS_BENCHMARK) # # Benchmarks # - - message ("searching: ${DEPS_ROOT} for google benchmark libs") - find_library (BENCHMARK_LIBRARY benchmark PATHS ${DEPS_ROOT}/lib) - find_library (BENCHMARKM_LIBRARY benchmark_main PATHS ${DEPS_ROOT}/lib) - #find_library (PTHREAD_LIBRARY pthread) - message ("BENCHMARK_LIBRARY: ${BENCHMARK_LIBRARY}") - include_directories ("${BENCHMARK_INC_DIR}") if (Boost_FOUND AND NO) add_definitions (-DHAVE_BOOST=1) include_directories ("${Boost_INCLUDE_DIRS}") @@ -110,11 +113,14 @@ if (CPPDUALS_BENCHMARK) set (BLA_STATIC OFF) endif (NOT BLA_STATIC) set (BLA_VENDOR OpenBLAS) - endif (NOT APPLE AND NOT BLA_VENDOR) + endif () find_package (BLAS REQUIRED) #find_package (LAPACK REQUIRED) add_definitions (-DHAVE_BLAS) #add_definitions (-DEIGEN_USE_BLAS) + if (APPLE) + add_definitions (-DACCELERATE_NEW_LAPACK -DACCELERATE_LAPACK_ILP64) + endif () # find lapacke.h cblas.h set (CBLAS_HINTS ${BLAS_DIR} ${LAPACK_DIR} /usr /usr/local /opt /opt/local) @@ -124,18 +130,21 @@ if (CPPDUALS_BENCHMARK) /opt /opt/local /usr/local/opt - /System/Library/Frameworks) + /System/Library/Frameworks + ${BLAS_LIBRARIES} + ) # Finds the include directories for lapacke.h find_path (LAPACKE_INCLUDE_DIRS - NAMES lapacke.h + REQUIRED + NAMES lapacke.h clapack.h HINTS ${CBLAS_HINTS} PATH_SUFFIXES include inc include/x86_64 include/x64 - openblas/include -# Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Versions/Current/Headers + openblas/include openblas + Frameworks/vecLib.framework/Headers PATHS ${CBLAS_PATHS} - DOC "LAPACK(E) include header lapacke.h") + DOC "LAPACK(E) include header lapacke.h/clapack.h") mark_as_advanced (LAPACKE_INCLUDE_DIRS) if (LAPACKE_INCLUDE_DIRS) include_directories (${LAPACKE_INCLUDE_DIRS}) @@ -145,12 +154,13 @@ if (CPPDUALS_BENCHMARK) # Finds the include directories for cblas*.h find_path (CBLAS_INCLUDE_DIRS + REQUIRED NAMES cblas.h cblas_openblas.h cblas-openblas.h HINTS ${CBLAS_HINTS} PATH_SUFFIXES include inc include/x86_64 include/x64 - openblas/include -# Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Versions/Current/Headers + openblas/include openblas + Frameworks/vecLib.framework/Headers PATHS ${CBLAS_PATHS} DOC "BLAS include header cblas.h") mark_as_advanced (CBLAS_INCLUDE_DIRS) @@ -161,80 +171,168 @@ if (CPPDUALS_BENCHMARK) break() endif (EXISTS "${CBLAS_INCLUDE_DIRS}/${cblas}") endforeach (cblas) - message ("Found BLAS : ${BLAS_LIBRARIES}") message ("Found cBLAS : ${CBLAS_INCLUDE_DIRS}") message ("Found lapacke : ${LAPACKE_INCLUDE_DIRS}") - set (OPT_FLAGS "") + set (BMK_FLAGS "") if (NOT MSVC) - #set (OPT_FLAGS "-O3;-mavx") - #set (OPT_FLAGS "-O3;-march=native;-fopenmp") - set (OPT_FLAGS "-O3;-msse3;-fopenmp") - set (OPT_FLAGS "-O3") - set (OPT_FLAGS "-O3;-msse3") - set (OPT_FLAGS "-O3;-march=native;-funroll-loops") - set (OPT_FLAGS "-O3;-msse3;-mavx2;-mfma") - set (OPT_FLAGS "-O3;-march=native") - #set (OPT_FLAGS "${OPT_FLAGS};-save-temps;-fverbose-asm") + #set (BMK_FLAGS "-O3;-mavx") + #set (BMK_FLAGS "-O3;-march=native;-fopenmp") + if (X86) + set (BMK_FLAGS "-msse3;-fopenmp") + set (BMK_FLAGS "-msse3") + set (BMK_FLAGS "-march=native;-funroll-loops") + set (BMK_FLAGS "-msse3;-mavx2;-mfma") + else () + set (BMK_FLAGS "-march=native") + endif () + #set (BMK_FLAGS "-O3;${BMK_FLAGS}") + #set (BMK_FLAGS "${BMK_FLAGS};-save-temps;-fverbose-asm") else () - set (OPT_FLAGS "/arch:IA32") - set (OPT_FLAGS "/arch:SSE") - set (OPT_FLAGS "/arch:SSE2") - set (OPT_FLAGS "/arch:AVX2") + set (BMK_FLAGS "/arch:IA32") + set (BMK_FLAGS "/arch:SSE") + set (BMK_FLAGS "/arch:SSE2") + set (BMK_FLAGS "/arch:AVX2") endif () - #set (OPT_FLAGS "${OPT_FLAGS};-DEIGEN_DONT_VECTORIZE") - #set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_DONT_VECTORIZE") - #set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_DONT_VECTORIZE_CDUAL") + #set (BMK_FLAGS "${BMK_FLAGS};-DEIGEN_DONT_VECTORIZE") + #set (BMK_FLAGS "${BMK_FLAGS};-DCPPDUALS_DONT_VECTORIZE") + #set (BMK_FLAGS "${BMK_FLAGS};-DCPPDUALS_DONT_VECTORIZE_CDUAL") - foreach (BENCH bench_dual bench_eigen bench_gemm bench_example bench_fmt) - add_executable (${BENCH} ${BENCH}.cpp) - target_compile_options (${BENCH} PUBLIC ${OPT_FLAGS}) - #set_target_properties (${BENCH} PROPERTIES LINK_FLAGS -fopenmp) - #target_link_options (${BENCH} PUBLIC ${OPT_FLAGS}) - string (REPLACE ";" ", " L2 "${OPT_FLAGS} ${CMAKE_CXX_FLAGS}") - target_compile_definitions (${BENCH} PRIVATE "OPT_FLAGS=${L2}") - add_dependencies (${BENCH} benchmarkX eigenX expokitX) - target_link_libraries (${BENCH} ${BENCHMARK_LIBRARY} -lpthread ${BLAS_LIBRARIES}) - endforeach (BENCH) + foreach (VECTORIZE YES NO) + foreach (BENCH bench_dual bench_eigen bench_exp bench_gemm bench_example bench_fmt) + if (NOT VECTORIZE) + set (BENCHE ${BENCH}) + else () + set (BENCHE ${BENCH}_novec) + endif () + add_executable (${BENCHE} ${BENCH}.cpp) + target_compile_options (${BENCHE} PUBLIC ${BMK_FLAGS}) + if (NOT VECTORIZE) + target_compile_options (${BENCHE} PUBLIC "-DEIGEN_DONT_VECTORIZE") + endif() + #set_target_properties (${BENCH} PROPERTIES LINK_FLAGS -fopenmp) + #target_link_options (${BENCH} PUBLIC ${BMK_FLAGS}) + string (REPLACE ";" ", " L2 "${BMK_FLAGS} ${CMAKE_CXX_FLAGS}") + target_compile_definitions (${BENCHE} PRIVATE "BMK_FLAGS=${L2}") + target_link_libraries (${BENCHE} benchmark::benchmark ${BLAS_LIBRARIES} eigen expokit) + endforeach () + endforeach () target_link_libraries (bench_fmt fmt::fmt) + target_link_libraries (bench_fmt_novec fmt::fmt) endif (CPPDUALS_BENCHMARK) add_executable (sandbox sandbox.cpp) -add_dependencies (sandbox eigenX expokitX ) # mpfrX mprealX -#target_compile_options (sandbox PUBLIC ${OPT_FLAGS}) +#target_compile_options (sandbox PUBLIC ${BMK_FLAGS}) target_compile_options (sandbox PUBLIC -DCPPDUALS_VECTORIZE_CDUAL) if (MSVC) - target_compile_options (sandbox PUBLIC /arch:AVX2) + if (X86) + target_compile_options (sandbox PUBLIC /arch:AVX2) + endif () else () - target_compile_options (sandbox PUBLIC -O1 -msse3 -mavx2 -mfma) + if (X86) + target_compile_options (sandbox PUBLIC -O1 -msse3 -mavx2 -mfma) + else () + target_compile_options (sandbox PUBLIC -O1 ) # -mfpu=neon + endif () endif () set_target_properties (sandbox PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}) -target_link_libraries (sandbox PUBLIC cppduals_coverage_config) +target_link_libraries (sandbox PUBLIC cppduals_coverage_config eigen expokit) # # Generate coverage reports # -if (CODE_COVERAGE) - add_custom_target (cov - DEPENDS ${ALL_TEST_BINS} - COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target test - COMMAND $ - COMMAND lcov --capture --directory . --output-file coverage.info - COMMAND lcov --remove coverage.info '/usr/*' --output-file coverage.info - COMMAND lcov --remove coverage.info '*/thirdparty/*' --output-file coverage.info - COMMAND lcov --remove coverage.info '*/googletest/*' --output-file coverage.info - COMMAND lcov --list coverage.info +if (CPPDUALS_CODE_COVERAGE) + # Find required tools + get_filename_component(CMAKE_CXX_COMPILER_DIR "${CMAKE_CXX_COMPILER}" DIRECTORY) + set(GENHTML_PATH ${CMAKE_CXX_COMPILER_DIR}/genhtml) + if(NOT EXISTS ${GENHTML_PATH}) + find_program(GENHTML_PATH genhtml REQUIRED) + endif() + + if (CMAKE_CXX_COMPILER_ID MATCHES "Clang") + # Clang + # often in the same directory as clang - search there too + get_filename_component(CMAKE_CXX_COMPILER_DIR "${CMAKE_CXX_COMPILER}" DIRECTORY) + set(LLVM_COV_PATH ${CMAKE_CXX_COMPILER_DIR}/llvm-cov) + set(LLVM_PROFDATA_PATH ${CMAKE_CXX_COMPILER_DIR}/llvm-profdata) + if(NOT EXISTS ${LLVM_COV_PATH}) + find_program(LLVM_COV_PATH llvm-cov REQUIRED) + endif() + if(NOT EXISTS ${LLVM_PROFDATA_PATH}) + find_program(LLVM_PROFDATA_PATH llvm-profdata REQUIRED) + endif() + + add_custom_target(cov + DEPENDS ${ALL_TEST_BINS} + COMMAND ${CMAKE_COMMAND} -E remove_directory coverage/profraw + COMMAND ${CMAKE_COMMAND} -E make_directory coverage/profraw + COMMAND ${CMAKE_COMMAND} --build . --target test "LLVM_PROFILE_FILE=coverage/profraw/%p.profraw" + COMMAND ${LLVM_PROFDATA_PATH} merge -sparse coverage/profraw/*.profraw -o coverage/coverage.profdata + COMMAND ${LLVM_COV_PATH} show + -instr-profile=coverage/coverage.profdata + -format=html -output-dir=coverage/html + -show-line-counts-or-regions + -show-instantiation-summary + ${ALL_TEST_BINS} + COMMAND ${LLVM_COV_PATH} report + -instr-profile=coverage/coverage.profdata + ${ALL_TEST_BINS} + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} ) - add_custom_target (cov-html + elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU") + # GCC + set(LCOV_PATH ${CMAKE_CXX_COMPILER_DIR}/lcov) + if(NOT EXISTS ${LCOV_PATH}) + find_program(LCOV_PATH lcov REQUIRED) + endif() + # Extract GCC version number and use matching gcov version + execute_process( + COMMAND ${CMAKE_CXX_COMPILER} -dumpversion + OUTPUT_VARIABLE GCC_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + set(GCOV_PATH ${CMAKE_CXX_COMPILER_DIR}/gcov-${GCC_VERSION}) + if(NOT EXISTS ${GCOV_PATH}) + find_program(GCOV_PATH gcov-${GCC_VERSION} REQUIRED) + endif() + + add_custom_target(cov + DEPENDS ${ALL_TEST_BINS} + COMMAND ${CMAKE_COMMAND} -E make_directory coverage + COMMAND ${LCOV_PATH} --directory . --zerocounters + # Add initial capture before running tests + COMMAND ${LCOV_PATH} --directory . --capture --initial --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.base --ignore-errors inconsistent + # Run tests and capture coverage data + COMMAND ctest --output-on-failure + COMMAND ${LCOV_PATH} --directory . --capture --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info --ignore-errors inconsistent + # Combine baseline and test coverage data + COMMAND ${LCOV_PATH} --add-tracefile coverage/coverage.base --add-tracefile coverage/coverage.info --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.total + # Only look at the coverage of the tests and duals library + COMMAND ${LCOV_PATH} --extract coverage/coverage.total '*/tests/*' '*/duals/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info + # Remove unwanted paths + #COMMAND ${LCOV_PATH} --remove coverage/coverage.total '/usr/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info + #COMMAND ${LCOV_PATH} --remove coverage/coverage.info '*/thirdparty/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info + #COMMAND ${LCOV_PATH} --remove coverage/coverage.info '*/googletest/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info + # Add --ignore-errors empty to prevent failure on empty coverage data + COMMAND ${LCOV_PATH} --list coverage/coverage.info --ignore-errors empty --gcov-tool ${GCOV_PATH} + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + ) + else() + message(FATAL_ERROR "No coverage tool found for ${CMAKE_CXX_COMPILER_ID}") + endif() + + add_custom_target(cov-html DEPENDS cov - COMMAND rm -rf ../coverage - COMMAND genhtml --ignore-errors source coverage.info --legend --title "make cov" - --output-directory=../coverage - COMMAND echo "output in coverage/index.html" - ) + COMMAND ${CMAKE_COMMAND} -E remove_directory coverage/html + COMMAND ${GENHTML_PATH} --ignore-errors source coverage/coverage.info --legend --title "cppduals coverage" + --output-directory=coverage/html + COMMAND ${CMAKE_COMMAND} -E echo "Coverage report generated at coverage/html/index.html" + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + ) + endif () diff --git a/src/include/cppduals/tests/bench_eigen.cpp b/src/include/cppduals/tests/bench_eigen.cpp index adde638d8..f32304276 100644 --- a/src/include/cppduals/tests/bench_eigen.cpp +++ b/src/include/cppduals/tests/bench_eigen.cpp @@ -11,13 +11,13 @@ // (c)2019 Michael Tesch. tesch1@gmail.com // +#include #include #include #include #include #include "type_name.hpp" -#include #include #include @@ -31,67 +31,6 @@ using namespace duals; template< class T > struct type_identity { typedef T type; }; -namespace Eigen { -namespace internal { -template struct is_exp_known_type; -template struct is_exp_known_type> : is_exp_known_type {}; -#if 0 -template struct MatrixExponentialScalingOp; -template -struct MatrixExponentialScalingOp> -{ - MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { } - inline const duals::dual operator() (const duals::dual & x) const - { - using std::ldexp; - return ldexp(x, -m_squarings); - } - typedef std::complex> ComplexScalar; - inline const ComplexScalar operator() (const ComplexScalar& x) const - { - using std::ldexp; - return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings)); - } - - private: - int m_squarings; -}; -#endif -}} -#include - -namespace Eigen { -namespace internal { -template -struct matrix_exp_computeUV > -{ - typedef typename NumTraits::Scalar>::Real RealScalar; - template - static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) - { - using std::frexp; - using std::pow; - const RealScalar l1norm = arg.cwiseAbs().colwise().sum().maxCoeff(); - squarings = 0; - if (l1norm < 1.495585217958292e-002) { - matrix_exp_pade3(arg, U, V); - } else if (l1norm < 2.539398330063230e-001) { - matrix_exp_pade5(arg, U, V); - } else if (l1norm < 9.504178996162932e-001) { - matrix_exp_pade7(arg, U, V); - } else if (l1norm < 2.097847961257068e+000) { - matrix_exp_pade9(arg, U, V); - } else { - const RealScalar maxnorm = 5.371920351148152; - frexp(l1norm / maxnorm, &squarings); - if (squarings < 0) squarings = 0; - MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp(squarings)); - matrix_exp_pade13(A, U, V); - } - } -}; -}} - /* encode the type into an integer for benchmark output */ template struct type_num { /* should fail */ }; template<> struct type_num { static constexpr int id = 1; }; @@ -271,87 +210,6 @@ template void B_MatVec(benchmark::State& state) { state.SetComplexityN(state.range(0)); } -template void B_Expm(benchmark::State& state) -{ - int N = state.range(0); - //Rt S(1); - MatrixX A = MatrixX::Random(N, N); - MatrixX B = MatrixX::Zero(N, N); - //A = S * A / A.norm(); - - for (auto _ : state) { - B = A.exp(); - benchmark::ClobberMemory(); - } - - state.counters["type"] = type_num::id; - state.SetComplexityN(state.range(0)); -} - -template void B_ExpPadm(benchmark::State& state) -{ - int N = state.range(0); - //Rt S(1); - MatrixX A = MatrixX::Random(N, N); - MatrixX B = MatrixX::Zero(N, N); - //A = S * A / A.norm(); - - for (auto _ : state) { - B = eexpokit::padm(A); - benchmark::ClobberMemory(); - } - - state.counters["type"] = type_num::id; - state.SetComplexityN(state.range(0)); -} - -template void B_ExpExpv(benchmark::State& state) -{ - int N = state.range(0); - //Rt S(1); - MatrixX A = MatrixX::Zero(N, N); - MatrixX b = MatrixX::Ones(N, 1); - MatrixX c = MatrixX::Zero(N, 1); - //A = S * A / A.norm(); - - // sparse random fill - for (int i = 0; i < 4*N; i++) - A((int)duals::randos::random(0.,N-1.), - (int)duals::randos::random(0.,N-1.)) = duals::randos::random2(); - - for (auto _ : state) { - auto ret = eexpokit::expv(1,A,b); - if (ret.err > 1) { - std::ofstream f("fail.m"); - f << "A=" << A.format(eexpokit::OctaveFmt) << "\n"; - break; - } - // c = ret.w - benchmark::ClobberMemory(); - } - - state.counters["type"] = type_num::id; - state.SetComplexityN(state.range(0)); -} - -template void B_ExpChbv(benchmark::State& state) -{ - int N = state.range(0); - //Rt S(1); - MatrixX A = MatrixX::Random(N, N); - MatrixX b = MatrixX::Zero(N, 1); - MatrixX c = MatrixX::Zero(N, 1); - //A = S * A / A.norm(); - - for (auto _ : state) { - c = eexpokit::chbv(A,b); - benchmark::ClobberMemory(); - } - - state.counters["type"] = type_num::id; - state.SetComplexityN(state.range(0)); -} - #define MAKE_BM_SIMPLE(TYPE1,TYPE2,NF) \ BENCHMARK_TEMPLATE(B_VecVecAdd, TYPE1,TYPE2) V_RANGE(4,NF); \ BENCHMARK_TEMPLATE(B_VecVecSub, TYPE1,TYPE2) V_RANGE(4,NF); \ @@ -359,11 +217,7 @@ template void B_ExpChbv(benchmark::State& state) BENCHMARK_TEMPLATE(B_VecVecDiv, TYPE1,TYPE2) V_RANGE(4,NF); \ BENCHMARK_TEMPLATE(B_MatVec, TYPE1,TYPE2) V_RANGE(4,NF); \ BENCHMARK_TEMPLATE(B_MatMat, TYPE1,TYPE2) V_RANGE(1,NF); \ - BENCHMARK_TEMPLATE(B_MatDiv, TYPE1,TYPE2) V_RANGE(1,NF); \ - BENCHMARK_TEMPLATE(B_Expm, TYPE1) V_RANGE(1,NF); \ - BENCHMARK_TEMPLATE(B_ExpPadm, TYPE1) V_RANGE(1,NF); \ - BENCHMARK_TEMPLATE(B_ExpChbv, TYPE1) V_RANGE(1,NF); \ - BENCHMARK_TEMPLATE(B_ExpExpv, TYPE1) V_RANGE(1,NF) + BENCHMARK_TEMPLATE(B_MatDiv, TYPE1,TYPE2) V_RANGE(1,NF) #define MAKE_BENCHMARKS(TYPE1,TYPE2,NF) \ MAKE_BM_SIMPLE(TYPE1,TYPE2,NF) @@ -389,7 +243,7 @@ MAKE_BM_SIMPLE(cduald, cduald,4); int main(int argc, char** argv) { #ifndef EIGEN_VECTORIZE - static_assert(false, "no vectorization?"); + //static_assert(false, "no vectorization?"); #endif std::cout << "OPT_FLAGS=" << QUOTE(OPT_FLAGS) << "\n"; std::cout << "INSTRUCTIONSET=" << Eigen::SimdInstructionSetsInUse() << "\n"; diff --git a/src/include/cppduals/tests/bench_gemm.cpp b/src/include/cppduals/tests/bench_gemm.cpp index 3426ccecc..53ed76b11 100644 --- a/src/include/cppduals/tests/bench_gemm.cpp +++ b/src/include/cppduals/tests/bench_gemm.cpp @@ -6,12 +6,14 @@ // // (c)2019 Michael Tesch. tesch1@gmail.com // +#include + #if defined(__APPLE__) && defined(__clang__) #include #else -#ifdef EIGEN_LAPACKE +#if defined(EIGEN_LAPACKE) || defined(__APPLE__) #include #else #include @@ -24,13 +26,13 @@ extern "C" { } #endif // defined(__APPLE__) && defined(__clang__) +#include #include #include #include #include #include "type_name.hpp" -#include #include #include @@ -80,6 +82,59 @@ template void B_MatMat(benchmark::State& state) { state.SetComplexityN(state.range(0)); } +// Helper templates to select correct BLAS routine +template +struct blas_gemm; + +#define TRANSPOSE(X) ((X) ? CblasTrans : CblasNoTrans) + +template <> +struct blas_gemm { + static void call(bool tA, bool tB, int m, int n, int k, float alpha, float *A, + int lda, float *B, int ldb, float beta, float *C, int ldc) + { + cblas_sgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, alpha, A, + lda, B, ldb, beta, C, ldc); + } +}; + +template <> +struct blas_gemm { + static void call(bool tA, bool tB, int m, int n, int k, double alpha, + double *A, int lda, double *B, int ldb, double beta, + double *C, int ldc) + { + cblas_dgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, alpha, A, + lda, B, ldb, beta, C, ldc); + } +}; + +template <> +struct blas_gemm> { + static void call(bool tA, bool tB, int m, int n, int k, + std::complex alpha, std::complex *A, int lda, + std::complex *B, int ldb, std::complex beta, + std::complex *C, int ldc) + { + cblas_cgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, &alpha, A, + lda, B, ldb, &beta, C, ldc); + } +}; + +template <> +struct blas_gemm> { + static void call(bool tA, bool tB, int m, int n, int k, + std::complex alpha, std::complex *A, int lda, + std::complex *B, int ldb, std::complex beta, + std::complex *C, int ldc) + { + cblas_zgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, &alpha, A, + lda, B, ldb, &beta, C, ldc); + } +}; + +#undef TRANSPOSE + template ::value>::type* = nullptr> void matrix_multiplcation(T *A, int Awidth, int Aheight, T *B, int Bwidth, int Bheight, @@ -99,31 +154,12 @@ void matrix_multiplcation(T *A, int Awidth, int Aheight, assert(A_width == B_height); int lda = tA ? m : k; int ldb = tB ? k : n; -#define TRANSPOSE(X) ((X) ? CblasTrans : CblasNoTrans) - // http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html - if (!is_complex::value) { - if (sizeof(T) == sizeof(float)) - cblas_sgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), - m, n, k, 1.0, (float *)A, lda, (float *)B, ldb, - std::real(beta), (float *)AB, n); - else - cblas_dgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), - m, n, k, 1.0, (double *)A, lda, (double *)B, ldb, - std::real(beta), (double *)AB, n); - } - else { - std::complex alphaf(1,0); - std::complex alpha(1,0); - if (Eigen::NumTraits::digits10() < 10) - cblas_cgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), - m, n, k, &alphaf, A, lda, B, ldb, &beta, AB, n); - else - cblas_zgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), - m, n, k, &alpha, A, lda, B, ldb, &beta, AB, n); - } -#undef TRANSPOSE + + // Call the appropriate BLAS routine based on type T + blas_gemm::call(tA, tB, m, n, k, T(1), A, lda, B, ldb, beta, AB, n); } + template ::value>::type* = nullptr> void matrix_multiplcation(T *A, int Awidth, int Aheight, T *B, int Bwidth, int Bheight, @@ -232,10 +268,10 @@ MAKE_BM_SIMPLE(cduald, cduald,4); int main(int argc, char** argv) { #ifndef EIGEN_VECTORIZE - static_assert(false, "no vectorization?"); + //static_assert(false, "no vectorization?"); #endif #ifndef NDEBUG - static_assert(false, "NDEBUG to benchmark?"); + //static_assert(false, "NDEBUG to benchmark?"); #endif std::cout << "OPT_FLAGS=" << QUOTE(OPT_FLAGS) << "\n"; std::cout << "INSTRUCTIONSET=" << Eigen::SimdInstructionSetsInUse() << "\n"; diff --git a/src/include/cppduals/tests/sandbox.cpp b/src/include/cppduals/tests/sandbox.cpp index e6b49b0e1..abae1e8d8 100644 --- a/src/include/cppduals/tests/sandbox.cpp +++ b/src/include/cppduals/tests/sandbox.cpp @@ -12,12 +12,12 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ +#include #include #include #include #include "type_name.hpp" -#include #include #include #include @@ -60,13 +60,45 @@ template struct common_type,T> { typedef Rando type; }; } -#if 0 +#if EIGEN_ARCH_ARM64 int main(int argc, char * argv[]) { emtx A,B,C; //emtx A,B,C; C = A * B; + + // + using namespace Eigen::internal; + float32x4_t a, b, v3; + float32x2_t f2; + float32x4x2_t yy; + float ff[] = {1., 2., 3., 4.}; + duals::dual x; + const uint64x2_t maskq = {0, (uint64_t)-1}; + const uint64x2_t imaskq = {(uint64_t)-1, 0}; + + std::cout << "a:" << std::hex << maskq[0] << "\n"; + std::cout << "b:" << maskq[1] << "\n"; + + std::cout << "v3=" << v3[0] << ", " << v3[1] << ", " << v3[2] << ", " << v3[3] << "\n"; + std::cout << "f2=" << f2[0] << ", " << f2[1] << "\n"; + + std::cout << ff[0] << ", " << ff[1] << ", " << ff[2] << ", " << ff[3] << "\n"; +#ifdef __ARM_NEON + std::cout << "__ARM_NEON:" << __ARM_NEON << "\n"; +#endif +#ifdef __ARM_NEON_FP + std::cout << "__ARM_NEON_FP:" << __ARM_NEON_FP << "\n"; +#endif +#ifdef __ARM_ARCH + std::cout << "__ARM_ARCH:" << __ARM_ARCH << "\n"; +#endif +#ifdef __ARM_BIG_ENDIAN + std::cout << "__ARM_BIG_ENDIAN:" << __ARM_BIG_ENDIAN << "\n"; +#endif + + return 0; } #elif 0 diff --git a/src/include/cppduals/tests/test_1.cpp b/src/include/cppduals/tests/test_1.cpp index 1203a2721..3f4bd1ce2 100644 --- a/src/include/cppduals/tests/test_1.cpp +++ b/src/include/cppduals/tests/test_1.cpp @@ -81,9 +81,7 @@ expm4(const Eigen::EigenBase & A_, R.setIdentity(); R += B; ReturnT S = B; - int ni = 0; for (int ii = 2; ii < maxt; ii++) { - ni++; S = Real(1.0/ii) * S * B; R += S; auto Sn = S.norm(); diff --git a/src/include/cppduals/tests/test_dual.cpp b/src/include/cppduals/tests/test_dual.cpp index 8c550b7d8..dbcc6c246 100644 --- a/src/include/cppduals/tests/test_dual.cpp +++ b/src/include/cppduals/tests/test_dual.cpp @@ -18,6 +18,7 @@ #include #include +#include #include "gtest/gtest.h" using duals::dualf; @@ -58,7 +59,7 @@ TEST(template_, dual_traits) { // depth EXPECT_EQ(dual_traits::depth, 0); EXPECT_EQ(dual_traits::depth, 0); - EXPECT_EQ(dual_traits::depth, 0); + EXPECT_EQ(dual_traits::depth, 1); EXPECT_EQ(dual_traits::depth, 1); EXPECT_EQ(dual_traits::depth, 2); } @@ -186,6 +187,8 @@ TEST(template_, common_type) { _EXPECT_TRUE(std::is_same); _EXPECT_TRUE(std::is_same); _EXPECT_TRUE(std::is_same); + _EXPECT_TRUE(std::is_same); + _EXPECT_TRUE(std::is_same); _EXPECT_TRUE(std::is_same); _EXPECT_TRUE(std::is_same); @@ -328,6 +331,7 @@ TEST(members, rpart) { EXPECT_EQ(z.rpart(), 4); EXPECT_EQ(z.dpart(), 3); } + TEST(members, dpart) { EXPECT_EQ(dpart(3), 0); dualf z(2,3); @@ -555,6 +559,15 @@ TEST(comparison, ge) { EXPECT_FALSE(1 >= a); } +#if 0 +TEST(simple_ops, add) { + // https://gitlab.com/tesch1/cppduals/-/issues/11 + duals::dual> x; + x = x+std::complex(1., 2.); + duals::dual> y = sqrt(x+std::complex(1., 2.)); +} +#endif + TEST(compound_assign, same_type) { // OP= dualf x = 2 + 4_e; @@ -977,25 +990,6 @@ TEST(non_class, random2) { EXPECT_NE(c2.dpart(), 0); EXPECT_NE(c1.rpart(), c2.rpart()); EXPECT_NE(c1.dpart(), c2.dpart()); - - // cdualf - cdualf d1 = duals::randos::random2(); - cdualf d2 = duals::randos::random2(); - EXPECT_NE(d1.real().rpart(), 0); - EXPECT_NE(d1.real().dpart(), 0); - EXPECT_NE(d1.imag().rpart(), 0); - EXPECT_NE(d1.imag().dpart(), 0); - - EXPECT_NE(d2.real().rpart(), 0); - EXPECT_NE(d2.real().dpart(), 0); - EXPECT_NE(d2.imag().rpart(), 0); - EXPECT_NE(d2.imag().dpart(), 0); - - EXPECT_NE(d1.real().rpart(), d2.real().rpart()); - EXPECT_NE(d1.real().dpart(), d2.real().dpart()); - - EXPECT_NE(d1.imag().rpart(), d2.imag().rpart()); - EXPECT_NE(d1.imag().dpart(), d2.imag().dpart()); } TEST(smoke, funcs) { @@ -1097,6 +1091,17 @@ TEST(complex, mixing) { // complex * real -> complex // complex * dual -> complex ? + // pow_dual_complex + B = pow(1_ef, C); + // pow_complex_dual + B = pow(C, 1_ef); + + // pow_dual_dual + B = pow(1_ef, 2_ef); + // pow_dual_scalar + B = pow(1_ef, 2); + // pow_scalar_dual + B = pow(2, 1_ef); } int main(int argc, char **argv) diff --git a/src/include/cppduals/tests/test_eigen.cpp b/src/include/cppduals/tests/test_eigen.cpp index 61dbaccd8..fa6cc5701 100644 --- a/src/include/cppduals/tests/test_eigen.cpp +++ b/src/include/cppduals/tests/test_eigen.cpp @@ -16,11 +16,12 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ -#include "type_name.hpp" #include +#include "type_name.hpp" #include #include #include + #include //#include #include "eexpokit/padm.hpp" @@ -39,10 +40,13 @@ using duals::is_complex; using duals::dual_traits; using namespace duals::literals; -typedef std::complex complexd; +typedef long double ldouble; typedef std::complex complexf; -typedef std::complex cduald; +typedef std::complex complexd; +typedef std::complex complexld; typedef std::complex cdualf; +typedef std::complex cduald; +typedef std::complex cdualld; template using emtx = Eigen::Matrix; template using smtx = Eigen::SparseMatrix; @@ -60,49 +64,6 @@ template using ecdf = Eigen::Matrix ; EXPECT_NEAR(rpart(A), rpart(B),tol); \ EXPECT_NEAR(dpart(A), dpart(B),tol) -/// Simple taylor series, truncated when |S| is "small enough" -template -ReturnT -expm4(const Eigen::EigenBase & A_, - typename DerivedA::RealScalar mn = std::numeric_limits::epsilon() * 1000) -{ - //std::cerr << "do PO:" << type_name() << "\n"; - typedef typename DerivedA::RealScalar Real; - using std::isfinite; - const DerivedA & A = A_.derived(); - int maxt = std::numeric_limits::digits; - int s = log2(rpart(A.derived().norm())) + 1; - s = std::max(0, s); - - auto B = A * pow(Real(2), -s); - ReturnT R(A.rows(), A.cols()); - R.setIdentity(); - R += B; - ReturnT S = B; - for (int ii = 2; ii < maxt; ii++) { - S = S * B * Real(1.0/ii); - R += S; - auto Sn = S.norm(); - if (!isfinite(Sn)) { - std::cout << "expm() non-finite norm:" << Sn << " at " << ii << "\n"; - std::cout << " |R| = " << R.norm() << " s=" << s << "\n" - << " |A| = " << rpart(A.real().norm()) << "\n" - << " |A/2^s|=" << rpart(A.real().norm()/pow(2,s)) << "\n"; - break; - } - // converged yet? - if (Sn < mn) - break; - if (ii == maxt - 1) { - std::cout << "expm() didn't converge in " << maxt << " |S| = " << Sn << "\n"; - throw std::invalid_argument("no converge"); - } - } - - for (; s; s--) - R = R * R; - return R; -} TEST(Eigen, NumTraits) { //::testing::StaticAssertTypeEq>::Real, float>(); @@ -448,11 +409,12 @@ TEST(measure, norm) { b = 3; Rt d(1); MatrixD x; - x << d; + x.array() = d; MatrixD a = (MatrixD() << 1,2,3, 4,5+5_ef,6, 7,8,9).finished(); //typename MatrixD::Index index; EXPECT_EQ(a.sum(), 45 + 5_ef); + EXPECT_EQ(x.sum(), 9); EXPECT_NEAR(rpart(a.norm()), 16.8819430161341337282, 1e-5); EXPECT_NEAR(rpart(a.mean()), 5, 1e-5); EXPECT_NEAR(dpart(a.mean()), 0.555555555555555, 1e-5); @@ -510,31 +472,49 @@ TEST(dpart, matrix) { EXPECT_EQ((dpart(AA) - CC).norm(),0); } -TEST(func, expm) { - typedef float T; - typedef dual dualt; - typedef std::complex> cdualt; -#define NN 3 -#define N2 6 - emtx a,b; - a = emtx::Random(); - a.array() += 1.1 + 2.2_ef; - a.setZero(); - a = eexpokit::padm(a); - EXPECT_LT((a - emtx::Identity()).norm(), 1e-6) << "a=" << a << "\n"; - a *= 1+2_e; - EXPECT_LT((a - emtx::Identity()).norm(), 1e-6) << "a=" << a << "\n"; +TEST(eigen, exp_typechecks) { + typedef emtx Mat; + typedef emtx Matd; + EXPECT_FALSE(Eigen::internal::is_exp_known_type::value); + EXPECT_FALSE(Eigen::internal::is_exp_known_type::value); - emtx c; - //b = a + 1_e * emtx::Random(); - c.setZero(); - c = c.exp(); - //c = expm4(c); - EXPECT_LT((c - emtx::Identity()).norm(), 1e-6) << "b=" << b << "\n"; - #undef NN - #undef N2 + typedef emtx Matc; + typedef emtx Matcd; + EXPECT_FALSE(Eigen::internal::is_exp_known_type::value); + EXPECT_FALSE(Eigen::internal::is_exp_known_type::value); } +const bool _exp = true; +const bool _padm = false; + +#define TEST_EXP(SCALAR_T, SIDE, EXP_OR_PADM) \ + TEST(func, exp_##SCALAR_T##_##SIDE##EXP_OR_PADM) { \ + emtx a,b; \ + a = emtx::Random(); \ + a.setZero(); \ + if (EXP_OR_PADM == _exp) a = a.exp(); \ + else a = eexpokit::padm(a); \ + EXPECT_LT((a - emtx::Identity()).norm(), 1e-6) << "a=" << a << "\n"; \ +} + +#define TEST_EXP_SIZES(SCALAR_T, EXP_OR_PADM) \ + TEST_EXP(SCALAR_T, 3, EXP_OR_PADM) \ + TEST_EXP(SCALAR_T, 4, EXP_OR_PADM) \ + TEST_EXP(SCALAR_T, 17, EXP_OR_PADM) + +#define TEST_EXP_SIZES_EOP(SCALAR_T) \ + TEST_EXP_SIZES(SCALAR_T, _padm) \ + TEST_EXP_SIZES(SCALAR_T, _exp) + +// just make sure padm is working +TEST_EXP_SIZES_EOP(float) +TEST_EXP_SIZES_EOP(double) +TEST_EXP_SIZES_EOP(ldouble) +TEST_EXP_SIZES_EOP(complexf) +TEST_EXP_SIZES_EOP(complexd) +TEST_EXP_SIZES_EOP(complexld) + +/* testing engine */ #define QUOTE(...) STRFY(__VA_ARGS__) #define STRFY(...) #__VA_ARGS__ diff --git a/src/include/cppduals/tests/test_expm.cpp b/src/include/cppduals/tests/test_expm.cpp index b42461d45..3e8e21e9b 100644 --- a/src/include/cppduals/tests/test_expm.cpp +++ b/src/include/cppduals/tests/test_expm.cpp @@ -16,8 +16,8 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ -#include "type_name.hpp" #include +#include "type_name.hpp" #include #include #include @@ -104,7 +104,7 @@ expm4(const Eigen::EigenBase & A_, return R; } -template > +template , int N2 = 2*NN> void dexpm() { //typedef std::complex T; //typedef std::complex> dualt; @@ -167,8 +167,6 @@ void dexpm() { << "eA2=" << eA2.block(0,0,std::min(4,NN),std::min(4,NN)) << "\n"; EXPECT_LT((dA1 - dA2).norm(), tol) << "dA1=" << dA1.block(0,0,std::min(4,NN),std::min(4,NN)) << "\n" << "dA2=" << dA2.block(0,0,std::min(4,NN),std::min(4,NN)) << "\n"; -#undef NN -#undef N2 } #if defined(PHASE_1) diff --git a/src/include/cppduals/tests/test_fmt.cpp b/src/include/cppduals/tests/test_fmt.cpp index 42111318d..376845e44 100644 --- a/src/include/cppduals/tests/test_fmt.cpp +++ b/src/include/cppduals/tests/test_fmt.cpp @@ -14,6 +14,7 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ #include +//#include #define CPPDUALS_LIBFMT #define CPPDUALS_LIBFMT_COMPLEX #include @@ -29,79 +30,131 @@ typedef std::complex cdualf; using namespace duals::literals; using namespace std::complex_literals; +TEST(libfmt, float_) { + std::string s = fmt::format("{:.1f}", 2.f); + EXPECT_EQ(s, "2.0"); +} +TEST(libfmt, double_) { + std::string s = fmt::format("{:.1f}", 2.); + EXPECT_EQ(s, "2.0"); +} TEST(libfmt, complex_) { - std::string s = fmt::format("{}", 2.f + 3if); + std::string s = fmt::format("{:.1f}", 2.f + 3if); EXPECT_EQ(s, "(2.0+3.0if)"); + s = fmt::format("{:.2f}", 2.f + 3if); + EXPECT_EQ(s, "(2.00+3.00if)"); } TEST(libfmt, complex_el) { std::string s = fmt::format("{}", 2.l + 3il); - EXPECT_EQ(s, "(2.0+3.0il)"); + EXPECT_EQ(s, "(2+3il)"); } TEST(libfmt, complex_flags) { std::string s; s = fmt::format("{}", 2. + 3i); - EXPECT_EQ(s, "(2.0+3.0i)"); - + EXPECT_EQ(s, "(2+3i)"); +} +TEST(libfmt, complex_flags_dollar) { + std::string s; s = fmt::format("{:$}", 2. + 3i); - EXPECT_EQ(s, "(2.0+3.0i)"); - + EXPECT_EQ(s, "(2+3i)"); +} +TEST(libfmt, complex_flags_star) { + std::string s; s = fmt::format("{:*}", 2. + 3i); + EXPECT_EQ(s, "(2+3*i)"); +} +TEST(libfmt, complex_flags_star_f) { + std::string s; + s = fmt::format("{:*.1f}", 2. + 3i); EXPECT_EQ(s, "(2.0+3.0*i)"); +} +TEST(libfmt, complex_flags_comma) { + std::string s; + s = fmt::format("{:,.1f}", 2. + 3i); + EXPECT_EQ(s, "(2.0,3.0)"); s = fmt::format("{:,}", 2. + 3i); - EXPECT_EQ(s, "(2.0,3.0)"); + EXPECT_EQ(s, "(2,3)"); +} + +TEST(libfmt, complex_flags_a) { + std::string s; // + + s = fmt::format("{:$+}", 2. + 3i); + EXPECT_EQ(s, "(+2+3i)"); + + s = fmt::format("{:$+.1f}", 2. + 3i); EXPECT_EQ(s, "(+2.0+3.0i)"); s = fmt::format("{:*+}", 2. + 3i); - EXPECT_EQ(s, "(+2.0+3.0*i)"); + EXPECT_EQ(s, "(+2+3*i)"); s = fmt::format("{:,+}", 2. + 3i); - EXPECT_EQ(s, "(+2.0,+3.0)"); + EXPECT_EQ(s, "(+2,+3)"); +} + +TEST(libfmt, complex_flags_b) { + std::string s; // + - s = fmt::format("{:$+}", 2. - 3i); - EXPECT_EQ(s, "(+2.0-3.0i)"); + EXPECT_EQ(s, "(+2-3i)"); s = fmt::format("{:*+}", 2. - 3i); - EXPECT_EQ(s, "(+2.0-3.0*i)"); + EXPECT_EQ(s, "(+2-3*i)"); s = fmt::format("{:,+}", 2. - 3i); + EXPECT_EQ(s, "(+2,-3)"); + + s = fmt::format("{:,+.1f}", 2. - 3i); EXPECT_EQ(s, "(+2.0,-3.0)"); } + TEST(libfmt, complex_all_real) { std::string s; s = fmt::format("{}", 2. + 0i); + EXPECT_EQ(s, "(2)"); + + s = fmt::format("{:.1f}", 2. + 0i); EXPECT_EQ(s, "(2.0)"); s = fmt::format("{:*}", 2. + 0i); - EXPECT_EQ(s, "(2.0)"); + EXPECT_EQ(s, "(2)"); s = fmt::format("{:,}", 2. + 0i); + EXPECT_EQ(s, "(2,0)"); + + s = fmt::format("{:,.1f}", 2. + 0i); EXPECT_EQ(s, "(2.0,0.0)"); } + TEST(libfmt, complex_all_imag) { std::string s; s = fmt::format("{}", 0. + 3i); + EXPECT_EQ(s, "(3i)"); + + s = fmt::format("{:.1f}", 0. + 3i); EXPECT_EQ(s, "(3.0i)"); s = fmt::format("{:*}", 0. + 3i); - EXPECT_EQ(s, "(3.0*i)"); + EXPECT_EQ(s, "(3*i)"); s = fmt::format("{:,}", 0. + 3i); - EXPECT_EQ(s, "(0.0,3.0)"); + EXPECT_EQ(s, "(0,3)"); } TEST(libfmt, complex_plus) { std::string s = fmt::format("{:+}", 1. + 3i); + EXPECT_EQ(s, "(+1+3i)"); + + s = fmt::format("{:+.1f}", 1. + 3i); EXPECT_EQ(s, "(+1.0+3.0i)"); s = fmt::format("{:+}", 1. - 3i); - EXPECT_EQ(s, "(+1.0-3.0i)"); + EXPECT_EQ(s, "(+1-3i)"); } TEST(libfmt, complex_g) { std::string s = fmt::format("{:g}", 2.f + 3if); @@ -118,67 +171,102 @@ TEST(libfmt, complex_gel) { TEST(libfmt, dual_) { std::string s = fmt::format("{}", 2 + 3_ef); + EXPECT_EQ(s, "(2+3_ef)"); + + s = fmt::format("{:.1f}", 2 + 3_ef); EXPECT_EQ(s, "(2.0+3.0_ef)"); } TEST(libfmt, dual_el) { std::string s = fmt::format("{}", 2 + 3_el); + EXPECT_EQ(s, "(2+3_el)"); + + s = fmt::format("{:.1f}", 2 + 3_el); EXPECT_EQ(s, "(2.0+3.0_el)"); } TEST(libfmt, dual_flags) { std::string s; s = fmt::format("{}", 2. + 3_e); - EXPECT_EQ(s, "(2.0+3.0_e)"); + EXPECT_EQ(s, "(2+3_e)"); s = fmt::format("{:$}", 2. + 3_e); + EXPECT_EQ(s, "(2+3_e)"); + + s = fmt::format("{:$.1f}", 2. + 3_e); EXPECT_EQ(s, "(2.0+3.0_e)"); s = fmt::format("{:*}", 2. + 3_e); - EXPECT_EQ(s, "(2.0+3.0*e)"); + EXPECT_EQ(s, "(2+3*e)"); s = fmt::format("{:,}", 2. + 3_e); + EXPECT_EQ(s, "(2,3)"); + + s = fmt::format("{:,.1f}", 2. + 3_e); EXPECT_EQ(s, "(2.0,3.0)"); // + + s = fmt::format("{:$+}", 2. + 3_e); - EXPECT_EQ(s, "(+2.0+3.0_e)"); + EXPECT_EQ(s, "(+2+3_e)"); s = fmt::format("{:*+}", 2. + 3_e); + EXPECT_EQ(s, "(+2+3*e)"); + + s = fmt::format("{:*+.1f}", 2. + 3_e); EXPECT_EQ(s, "(+2.0+3.0*e)"); s = fmt::format("{:,+}", 2. + 3_e); - EXPECT_EQ(s, "(+2.0,+3.0)"); + EXPECT_EQ(s, "(+2,+3)"); // + - s = fmt::format("{:$+}", 2. - 3_e); - EXPECT_EQ(s, "(+2.0-3.0_e)"); + EXPECT_EQ(s, "(+2-3_e)"); s = fmt::format("{:*+}", 2. - 3_e); + EXPECT_EQ(s, "(+2-3*e)"); + + s = fmt::format("{:*+.1f}", 2. - 3_e); EXPECT_EQ(s, "(+2.0-3.0*e)"); s = fmt::format("{:,+}", 2. - 3_e); - EXPECT_EQ(s, "(+2.0,-3.0)"); + EXPECT_EQ(s, "(+2,-3)"); + s = fmt::format("{:,+}", 2. + 3_e); + EXPECT_EQ(s, "(+2,+3)"); } TEST(libfmt, dual_all_real) { std::string s; s = fmt::format("{}", 2 + 0_e); - EXPECT_EQ(s, "(2.0)"); + EXPECT_EQ(s, "(2)"); s = fmt::format("{:*}", 2 + 0_e); - EXPECT_EQ(s, "(2.0)"); + EXPECT_EQ(s, "(2)"); s = fmt::format("{:,}", 2 + 0_e); - EXPECT_EQ(s, "(2.0,0.0)"); + EXPECT_EQ(s, "(2,0)"); } TEST(libfmt, dual_all_dual) { std::string s; s = fmt::format("a{}b", 0 + 3_e); + EXPECT_EQ(s, "a(3_e)b"); + s = fmt::format("a{:.1f}b", 0 + 3_e); EXPECT_EQ(s, "a(3.0_e)b"); s = fmt::format("a{:*}b", 0 + 3_e); - EXPECT_EQ(s, "a(3.0*e)b"); + EXPECT_EQ(s, "a(3*e)b"); s = fmt::format("a{:,}b", 0 + 3_e); + EXPECT_EQ(s, "a(0,3)b"); + s = fmt::format("a{:,.1f}b", 0 + 3_e); EXPECT_EQ(s, "a(0.0,3.0)b"); } + +TEST(libfmt, dual_plus) { + std::string s = fmt::format("{:+}", 1. + 3_e); + EXPECT_EQ(s, "(+1+3_e)"); + + s = fmt::format("{:+.1f}", 1. + 3_e); + EXPECT_EQ(s, "(+1.0+3.0_e)"); + + s = fmt::format("{:+}", 1. - 3_e); + EXPECT_EQ(s, "(+1-3_e)"); +} TEST(libfmt, dual_g) { std::string s = fmt::format("{:g}", 2 + 3_ef); EXPECT_EQ(s, "(2+3_ef)"); @@ -194,6 +282,12 @@ TEST(libfmt, dual_gel) { std::string s = fmt::format("{:g}", 2 + 3_el); EXPECT_EQ(s, "(2+3_el)"); } + +TEST(libfmt, dual_cd) { + std::string s = fmt::format("{}", cdualf(2 + 3_ef, 4 + 5_ef)); + EXPECT_EQ(s, "((2+3_ef)+(4+5_ef)i)"); +} + TEST(libfmt, dual_cgt) { std::string s = fmt::format("{:g}", cdualf(2 + 3_ef, 4 + 5_ef)); EXPECT_EQ(s, "((2+3_ef)+(4+5_ef)i)"); @@ -208,6 +302,10 @@ TEST(libfmt, dual_cgts) { s = fmt::format("{:,*g}", cdualf(2 + 3_ef, 4 + 5_ef)); EXPECT_EQ(s, "((2+3*ef),(4+5*ef))"); + // nonsense but should still work + s = fmt::format("{:*,g}", cdualf(2 + 3_ef, 4 + 5_ef)); + EXPECT_EQ(s, "((2,3)+(4,5)*i)"); + s = fmt::format("{:,,g}", cdualf(2 + 3_ef, 4 + 5_ef)); EXPECT_EQ(s, "((2,3),(4,5))"); } diff --git a/src/include/cppduals/tests/test_funcs.cpp b/src/include/cppduals/tests/test_funcs.cpp index c897cad45..ae9947214 100644 --- a/src/include/cppduals/tests/test_funcs.cpp +++ b/src/include/cppduals/tests/test_funcs.cpp @@ -32,25 +32,38 @@ typedef std::complex complexd; // not verify precision. #define FD_CHECK(T, F, ...) \ TEST(func##_##T, F) { \ + using std::isnan; \ for (T x : __VA_ARGS__) { \ T prec = 100 * std::sqrt(std::numeric_limits::epsilon()); \ - T dd = dpart(F(x + dual(0,1))); \ + T dd = duals::dpart(F(x + dual(0,1))); \ /*T dx = std::numeric_limits::epsilon() * (T)1000000; */ \ T dx = T(1)/ (1ull << (std::numeric_limits::digits / 3)); \ T fd = (F(x + dx) - F(x - dx)) / (2*dx); \ - EXPECT_CNEAR(dd, fd, prec * std::abs(std::max(std::max(dd,fd),T(1)))) \ - << "dd=" << dd << " fd=" << fd << " x=" << x << " dx=" << dx; \ + if (!isnan(dd) && !isnan(fd)) { \ + EXPECT_CNEAR(dd, fd, \ + prec * std::abs(std::max(std::max(dd,fd),T(1)))) \ + << "dd=" << dd << " fd=" << fd << " x=" << x << " dx=" << dx; \ + } \ } \ } +#define powL(x) pow(x,2) +#define powR(x) pow(2,x) +#define powLR(x) pow(x,x) + FD_CHECK(double, exp, {-1,0,1}) -FD_CHECK(double, log, {1}) -//FD_CHECK(complexd, log, {-1,1}) TODO -FD_CHECK(double, log10, {1}) -//FD_CHECK(complexd, log10, {-1,0,1}) TODO +FD_CHECK(double, log, {1, 3}) +// FD_CHECK(complexd, log, {-1,1}) TODO +FD_CHECK(double, log10, {1, 3}) + +FD_CHECK(double, powL, {-3.,-1.,-0.4,0.,0.6,1.,2.5}) +FD_CHECK(double, powR, {-3.,-1.,-0.4,0.,0.6,1.,2.5}) +FD_CHECK(double, powLR, {-3.,-1.,-0.4,0.,0.6,1.,2.5}) + +// FD_CHECK(complexd, log10, {-1,0,1}) TODO FD_CHECK(double, sqrt, {0.5,1.0}) FD_CHECK(double, cbrt, {-10.,-0.01,0.01,1.0,10.}) -//FD_CHECK(complexd, sqrt, {0,1}) TODO +// FD_CHECK(complexd, sqrt, {0,1}) TODO FD_CHECK(double, sin, {-1,0,1}) FD_CHECK(double, cos, {-1,0,1}) FD_CHECK(double, tan, {-1,0,1}) @@ -59,12 +72,27 @@ FD_CHECK(double, acos, {-.9,0.,.9}) FD_CHECK(double, atan, {-10,-1,0,1,10}) // TODO: -//FD_CHECK(double, sinh, {0}) -//FD_CHECK(double, cosh, {0}) -//FD_CHECK(double, tanh, {0}) -//FD_CHECK(double, asinh, {0}) -//FD_CHECK(double, acosh, {0}) -//FD_CHECK(double, atanh, {0}) +#define atan2L(x) atan2(x,2) +#define atan2R(x) atan2(2,x) +#define atan2LR(x) atan2(x,x) +FD_CHECK(double, atan2L, {-10.,-1.,0.,1.,10.}) +FD_CHECK(double, atan2R, {-10.,-1.,0.01,1.,10.}) +FD_CHECK(double, atan2LR, {-10.,-1.,0.01,1.,10.}) + +#define hypot2LR(x) hypot(x,x) +FD_CHECK(double, hypot2LR, {-10.,-1.,0.01,1.,10.}) + +#define scalbnL(x) scalbn(x,2) +FD_CHECK(double, scalbnL, {-10.,-1.,0.01,1.,10.}) + +FD_CHECK(double, sinh, {-0.1, 0.1}) +FD_CHECK(double, cosh, {-0.1, 0.1}) +FD_CHECK(double, tanh, {-0.1, 0.1}) +FD_CHECK(double, asinh, {-0.1, 0.1}) +FD_CHECK(double, acosh, {-1.1, 1.1}) +FD_CHECK(double, atanh, {-0.1, 0.1}) +FD_CHECK(double, log1p, {-0.1, 0.1}) +FD_CHECK(double, expm1, {-0.1, 0.1}) FD_CHECK(double, erf, {-1,0,1}) FD_CHECK(double, erfc, {-1,0,1}) @@ -97,6 +125,8 @@ TEST(func, tgamma) { //EXPECT_EQ(tgamma(x).rpart(), 362880); "interestingly", compiling without optimization (-O0) causes this to fail EXPECT_NEAR(tgamma(x).rpart(), 362880, 362880 * 100 * std::numeric_limits::epsilon()); } + +// part selection functions TEST(func, rpart) { dualf x = 10 + 4_e; EXPECT_EQ(rpart(x), 10); @@ -105,16 +135,258 @@ TEST(func, dpart) { dualf x = 2 + 4_e; EXPECT_EQ(dpart(x), 4); } + +// non-differentiable operations on the real part. TEST(func, abs) { + dualf x = -10 + 4_e; + EXPECT_EQ(abs(x), 10); } -TEST(func, arg) { +TEST(func, fabs) { + dualf x = -10 + 4_e; + EXPECT_EQ(fabs(x), 10); +} +TEST(func, fmax) { + dualf x = -10 + 4_e; + dualf y = 10 + 4_e; + EXPECT_EQ(fmax(x, y), 10); +} +TEST(func, fmin) { + dualf x = -10 + 4_e; + dualf y = 10 + 4_e; + EXPECT_EQ(fmin(x, y), -10); +} +TEST(func, frexp) { + dualf x = 6 + 4_e; + int exp = 0; + EXPECT_EQ(frexp(x, &exp), 0.75); + EXPECT_EQ(exp, 3); +} +TEST(func, ldexp) { + dualf x = 0.5 + 4_e; + int exp = 1; + EXPECT_EQ(ldexp(x, exp), 1); +} +TEST(func, trunc) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(trunc(x), 1); +} +TEST(func, floor) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(floor(x), 1); +} +TEST(func, ceil) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(ceil(x), 2); +} +TEST(func, round) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(round(x), 2); +} +// floating point functions +TEST(func, fpclassify) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(fpclassify(x), FP_NORMAL); + EXPECT_EQ(fpclassify(std::numeric_limits::infinity()), FP_INFINITE); + EXPECT_EQ(fpclassify(std::numeric_limits::quiet_NaN()), FP_NAN); + if (std::numeric_limits::has_denorm != std::denorm_absent) { + EXPECT_EQ(fpclassify(std::numeric_limits::denorm_min()), FP_SUBNORMAL); + } + EXPECT_EQ(fpclassify(x+std::numeric_limits::min()), FP_NORMAL); + EXPECT_EQ(fpclassify(2*std::numeric_limits::max()), FP_INFINITE); + EXPECT_EQ(fpclassify(x+std::numeric_limits::epsilon()), FP_NORMAL); +} +TEST(func, isfinite) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(isfinite(x), true); +} +TEST(func, isnormal) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(isnormal(x), true); +} +TEST(func, isinf) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(isinf(x), false); +} +TEST(func, isnan) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(isnan(x), false); +} +TEST(func, signbit) { + dualf x = 1.5 + 4_e; + EXPECT_EQ(signbit(x), false); +} +TEST(func, copysign) { + dualf x = 1.5 + 4_e; + dualf y = -1.3 + 2_e; + EXPECT_EQ(copysign(x, y), -1.5); +} + +// Utility functions +TEST(func, random) { + dualf x = random(0.001 + 0.001_e, 1 + 1_e); + EXPECT_GE(rpart(x), 0.001); + EXPECT_LE(rpart(x), 1); + EXPECT_GE(dpart(x), 0.001); + EXPECT_LE(dpart(x), 1); +} + +TEST(func, random2) { + dualf x = duals::randos::random2(0.001 + 0.001_e, 1 + 1_e); + EXPECT_GE(rpart(x), 0.001); + EXPECT_LE(rpart(x), 1); + EXPECT_GE(dpart(x), 0.001); + EXPECT_LE(dpart(x), 1); +} + +// more tests +TEST(func, logb) { + dualf x = 4 + 1_e; + EXPECT_EQ(rpart(logb(.5 + 1_e)), -1.); + EXPECT_EQ(rpart(logb(1 + 1_e)), 0.); + EXPECT_EQ(rpart(logb(2 + 1_e)), 1.); + EXPECT_EQ(rpart(logb(3 + 1_e)), 1.); + EXPECT_EQ(rpart(logb(4 + 1_e)), 2.); + EXPECT_EQ(rpart(logb(x * x)), 4.); + + EXPECT_EQ(dpart(logb(4 + 8_e)), std::numeric_limits::infinity()); + EXPECT_EQ(dpart(logb(4.01 + 8_e)), 0.); + + EXPECT_EQ(dpart(logb(x * x)), std::numeric_limits::infinity()); + EXPECT_EQ(dpart(logb(3 * x)), 0.); + x += 0.01; + EXPECT_EQ(dpart(logb(3 * x)), 0.); + EXPECT_EQ(dpart(logb(x * x)), 0.); +} + +TEST(func, pow) { + dualf x = pow(0 + 0_e, 0.); + EXPECT_EQ(rpart(x), 1); + EXPECT_EQ(dpart(x), 0); + + dualf y = pow(0 + 0_e, 0. + 0_e); + EXPECT_EQ(rpart(y), 1); + EXPECT_EQ(dpart(y), 0); + + dualf z = pow(0, 0. + 0_e); + EXPECT_EQ(rpart(z), 1); + EXPECT_EQ(dpart(z), 0); +} +TEST(func, pow_complex) { + std::complex C(3+4_ef, 5+6_ef); + std::complex x = std::pow(C, 2+1_ef); + std::complex ref = std::complex(std::pow(std::complex(3,5), 2)); + EXPECT_NEAR(rpart(x.real()), ref.real(), 1e-5); + EXPECT_NEAR(rpart(x.imag()), ref.imag(), 1e-6); + //EXPECT_EQ(dpart(x), std::pow(std::complex(3, 5), 2) * (2+1_ef)); +} + +//---------------------------------------------------------------------- +// Test for pow_dual_complex(const dual& realBase, const std::complex>& complexExponent) +//---------------------------------------------------------------------- +TEST(func, PowDualComplexTest) +{ + // 1) Prepare inputs: + // realBase = 2.0 (with zero derivative for this example) + duald realBase(2.0f); + + // complexExponent = (1.5 + 0.7i), + // each part a dual with .rpart() = 1.5 or 0.7, derivative 0 + std::complex complexExponent(duald(1.5f), duald(0.7f)); + + // 2) Call the function we want to test: + // x^y => pow_dual_complex(realBase, complexExponent) + std::complex result = std::pow(realBase, complexExponent); + + // 3) Reference: use standard pow in double + double dblBase = 2.0; + std::complex dblExponent(1.5, 0.7); + std::complex reference = std::pow(dblBase, dblExponent); + + // 4) Compare the .rpart() of the dual’s real/imag with reference + EXPECT_NEAR(result.real().rpart(), reference.real(), 1e-6); + EXPECT_NEAR(result.imag().rpart(), reference.imag(), 1e-6); +} + +//---------------------------------------------------------------------- +// Test for pow_complex_dual(const std::complex>& complexBase, const dual& realExponent) +//---------------------------------------------------------------------- +TEST(func, PowComplexDualTest) +{ + // 1) Prepare inputs: + // complexBase = (3 + 5i), each part dual with zero derivative + std::complex complexBase(duald(3.0f), duald(5.0f)); + + // realExponent = 2.0 as a dual + duald realExponent(2.0f); + + // 2) Call the function we want to test: + // x^y => pow_complex_dual(complexBase, realExponent) + std::complex result = pow(complexBase, realExponent); + + // 3) Reference: again, standard pow in double + std::complex dblBase(3.0, 5.0); + double dblExponent = 2.0; + std::complex reference = std::pow(dblBase, dblExponent); + + // 4) Compare the .rpart() of the dual’s real/imag with reference + EXPECT_NEAR(result.real().rpart(), reference.real(), 1e-6); + EXPECT_NEAR(result.imag().rpart(), reference.imag(), 1e-6); } TEST(func, norm) { + // TODO } TEST(func, conj) { + // TODO } TEST(func, polar) { + // TODO } +TEST(func, atan) { + EXPECT_EQ(rpart(atan(0 + 1_e)), atan(0)); + EXPECT_EQ(dpart(atan(1_e)), 1); + EXPECT_EQ(dpart(atan(1 + 1_e)), 0.5); // = 1 / (1 + x^2) + EXPECT_EQ(dpart(atan(-2 + 1_e)), 1. / 5.); // = 1 / (1 + x^2) +} +TEST(func, atan2) { + // TODO + //EXPECT_EQ(dpart(atan2(1_e, 1)), 1); + //EXPECT_EQ(dpart(atan2(1 + 1_e, 1)), 0.5); + //EXPECT_EQ(dpart(atan2(-2 + 1_e, 1)), (1. / 5.)); + duald y = 1 + 1_e; + duald x = 1 + 0_e; + auto z = atan2(y, x); + z = atan2(y, x); + EXPECT_EQ(rpart(z), atan2(rpart(y), rpart(x))); + EXPECT_EQ(dpart(z), 0.5); + + y = -2 + 1_e; + x = 1 + 0_e; + z = atan2(y, x); + EXPECT_EQ(rpart(z), atan2(rpart(y), rpart(x))); + EXPECT_EQ(dpart(z), 1. / 5.); + + y = 1 + 0_e; + x = -2 + 1_e; + z = atan2(y, x); + EXPECT_EQ(rpart(z), atan2(rpart(y), rpart(x))); + EXPECT_EQ(dpart(z), -1. / 5.); +} +TEST(func, atan2a) { + // TODO + duald y = 1 + 1_e; + EXPECT_EQ(rpart(atan2(y, 1)), atan2(rpart(y), 1)); + EXPECT_EQ(rpart(atan2(y, 2)), atan2(rpart(y), 2)); + EXPECT_EQ(dpart(atan2(y, 1)), 0.5); + y = 2 + 1_e; + EXPECT_EQ(dpart(atan2(y, -2)), -0.25); +} +TEST(func, atan2b) { + // TODO + duald x = 10 + 1_e; + EXPECT_EQ(rpart(atan2(2, x)), atan2(2, rpart(x))); + EXPECT_EQ(dpart(atan2(2, x)), -1./52); +} + struct pike_f1 { // function diff --git a/src/include/cppduals/tests/test_packets.cpp b/src/include/cppduals/tests/test_packets.cpp index c6f925a12..6300095ae 100644 --- a/src/include/cppduals/tests/test_packets.cpp +++ b/src/include/cppduals/tests/test_packets.cpp @@ -16,8 +16,8 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ -#include "type_name.hpp" #include +#include "type_name.hpp" #include #include #include @@ -60,6 +60,7 @@ template using ecdf = Eigen::Matrix ; #if !defined(CPPDUALS_DONT_VECTORIZE) && !defined(EIGEN_DONT_VECTORIZE) +#ifdef EIGEN_VECTORIZE_SSE TEST(Packet1cdf, pload_pstore) { using namespace Eigen::internal; cdualf cd1 = cdualf(1+2_ef,3+4_ef); @@ -68,6 +69,7 @@ TEST(Packet1cdf, pload_pstore) { pstore(&cd2, p1); EXPECT_DEQ(cd1, cd2); } +#endif using duals::randos::random2; @@ -89,11 +91,11 @@ using duals::randos::random2; pstore(cd3.data(), p3); \ for (int i = 0; i < N; i++) { \ EXPECT_DNEAR(cd3[i], cd1[i] op cd2[i], tol) \ - << cd1[i] << ',' << cd2[i] << " fail at " << i; \ + << cd1[i] << #op << cd2[i] << "=" << cd3[i] << "!=" << (cd1[i] op cd2[i]) << " fail at " << i << "/" << N; \ } \ } -#define GEN_PACKET_TEST_UN(PTYPE,pop,op) TEST(PTYPE,pop) { \ +#define GEN_PACKET_TEST_UN(PTYPE,pop,op) TEST(PTYPE,pop) { \ using namespace Eigen::internal; \ typedef unpacket_traits::type DTYPE; \ const static int N = unpacket_traits::size; \ @@ -112,7 +114,7 @@ using duals::randos::random2; } \ } -#define GEN_PACKET_TEST_RD(PTYPE,pop,op) TEST(PTYPE,pop) { \ +#define GEN_PACKET_TEST_RD(PTYPE,pop,op) TEST(PTYPE,pop) { \ using namespace Eigen::internal; \ typedef unpacket_traits::type DTYPE; \ const static int N = unpacket_traits::size; \ @@ -131,7 +133,7 @@ using duals::randos::random2; << acc << " " << p2 << " fail."; \ } -#define GEN_PACKET_TEST_REVERSE(PTYPE) TEST(PTYPE,preverse) { \ +#define GEN_PACKET_TEST_REVERSE(PTYPE) TEST(PTYPE,preverse) { \ using namespace Eigen::internal; \ typedef unpacket_traits::type DTYPE; \ const static int N = unpacket_traits::size; \ @@ -149,7 +151,7 @@ using duals::randos::random2; } \ } -#define GEN_PACKET_TEST_FIRST(PTYPE) TEST(PTYPE,pfirst) { \ +#define GEN_PACKET_TEST_FIRST(PTYPE) TEST(PTYPE,pfirst) { \ using namespace Eigen::internal; \ typedef unpacket_traits::type DTYPE; \ const static int N = unpacket_traits::size; \ @@ -323,13 +325,11 @@ using duals::randos::random2; GEN_PACKET_TEST_UN(PTYPE,pconj,conj) // TODO: -//pcplxflip //preduxp //pand //por //pxor //andnot -//pbroadcast4 //ploadquad (for packets w/ size==8) //pgather //pscatter @@ -345,7 +345,7 @@ GEN_CPACKET_TESTS(Packet2cf) GEN_CPACKET_TESTS(Packet4cf) #endif -#ifdef EIGEN_VECTORIZE_SSE +#if defined(EIGEN_VECTORIZE_SSE) || defined(EIGEN_VECTORIZE_NEON) GEN_PACKET_TESTS(Packet2df) GEN_PACKET_TESTS(Packet1dd) GEN_CPACKET_TESTS(Packet1cdf) @@ -368,20 +368,28 @@ TEST(compile, VECTORIZE) { #endif } TEST(compile, SSE) { +#ifdef EIGEN_VECTORIZE_SSE EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("SSE") != std::string::npos) - << "Not using SSE instructions:" << Eigen::SimdInstructionSetsInUse(); -#ifndef EIGEN_VECTORIZE_SSE - EXPECT_TRUE(false) - << "Not using EIGEN_VECTORIZE_SSE:" << Eigen::SimdInstructionSetsInUse(); +#else + EXPECT_TRUE(true) #endif + << "Not using EIGEN_VECTORIZE_SSE:" << Eigen::SimdInstructionSetsInUse(); } TEST(compile, AVX) { +#ifdef EIGEN_VECTORIZE_AVX EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("AVX") != std::string::npos) - << "Not using AVX instructions:" << Eigen::SimdInstructionSetsInUse(); -#ifndef EIGEN_VECTORIZE_AVX - EXPECT_TRUE(false) - << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse(); +#else + EXPECT_TRUE(true) #endif + << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse(); +} +TEST(compile, NEON) { +#ifdef EIGEN_VECTORIZE_NEON + EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("NEON") != std::string::npos) +#else + EXPECT_TRUE(true) +#endif + << "Not using EIGEN_VECTORIZE_NEON:" << Eigen::SimdInstructionSetsInUse(); } #define QUOTE(...) STRFY(__VA_ARGS__) diff --git a/src/include/cppduals/tests/test_solve.cpp b/src/include/cppduals/tests/test_solve.cpp index a7c9c8008..326a273f9 100644 --- a/src/include/cppduals/tests/test_solve.cpp +++ b/src/include/cppduals/tests/test_solve.cpp @@ -16,8 +16,8 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ -#include "type_name.hpp" #include +#include "type_name.hpp" #include #include #include @@ -78,7 +78,10 @@ void solveLu() { emtx B = b + DT(0,1) * emtx::Random(); emtx C,D,E; C = A * B; - D = A.lu().solve(C); + // D = A.lu().solve(C); + // D = A.colPivHouseholderQr().solve(C); + // D = A.bdcSvd().solve(C); + D = A.fullPivLu().solve(C); EXPECT_LT(rpart(B - D).norm(), tol); EXPECT_LT(dpart(B - D).norm(), tol); } diff --git a/src/include/cppduals/tests/test_vectorize.cpp b/src/include/cppduals/tests/test_vectorize.cpp index 1299c7895..0e3560a37 100644 --- a/src/include/cppduals/tests/test_vectorize.cpp +++ b/src/include/cppduals/tests/test_vectorize.cpp @@ -16,8 +16,8 @@ * (c)2019 Michael Tesch. tesch1@gmail.com */ -#include "type_name.hpp" #include +#include "type_name.hpp" #include #include #include @@ -52,7 +52,7 @@ template using ecdf = Eigen::Matrix ; #define _EXPECT_FALSE(...) {typedef __VA_ARGS__ fal; EXPECT_FALSE(fal::value); static_assert(!fal::value, "sa"); } #define ASSERT_DEQ(A,B) ASSERT_EQ(rpart(A), rpart(B)); ASSERT_EQ(dpart(A), dpart(B)) #define ASSERT_DNEAR(A,B,tol) \ - ASSERT_NEAR(abs(rpart((A) - (B))),0,abs(rpart(A))*(tol)); \ + ASSERT_NEAR(abs(rpart((A) - (B))),0,abs(rpart(A))*(tol)) << "rpart " << A << " " << B << "\n"; \ ASSERT_NEAR(abs(dpart((A) - (B))),0,abs(dpart(A))*(tol)) #define EXPECT_DEQ(A,B) EXPECT_EQ(rpart(A), rpart(B)); EXPECT_EQ(dpart(A), dpart(B)) #define EXPECT_DNE(A,B) EXPECT_NE(rpart(A), rpart(B)); EXPECT_NE(dpart(A), dpart(B)) @@ -117,9 +117,9 @@ void elemwise(int N) { ASSERT_DEQ(cca[i], Cca(i)) << "ca mismatch at " << i << "\n"; ASSERT_DEQ(ccb[i], Ccb(i)) << "cb mismatch at " << i << "\n"; } - ASSERT_DNEAR(sum, A.sum(), N*tol); - ASSERT_DNEAR(sum, A.sum(), N*tol); + ASSERT_DNEAR(sum, A.sum(), N*tol) << "sum mismatch at " << sum << " " << A.sum() << "\n"; } +#if defined(PHASE_1) TEST(Vector, full_even_dualf) { elemwise(512); } TEST(Vector, full_even_duald) { elemwise(512); } @@ -136,6 +136,7 @@ TEST(Vector, single_elem_cdualf) { elemwise(1); } TEST(Vector, two_elem_dualf) { elemwise(2); } TEST(Vector, two_elem_duald) { elemwise(2); } TEST(Vector, two_elem_cdualf) { elemwise(2); } +#endif #define DBOUT(X) #define MAKE_MULT_TEST(TYPE1, TYPE2, FIX, SIZE) \ @@ -308,20 +309,30 @@ TEST(flags, VECTORIZE) { #endif } TEST(flags, SSE) { +#ifdef EIGEN_VECTORIZE_SSE EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("SSE") != std::string::npos) - << "Not using SSE instructions:" << Eigen::SimdInstructionSetsInUse(); -#ifndef EIGEN_VECTORIZE_SSE - EXPECT_TRUE(false) +#else + EXPECT_TRUE(true) +#endif << "Not using EIGEN_VECTORIZE_SSE:" << Eigen::SimdInstructionSetsInUse(); -#endif } + TEST(flags, AVX) { +#ifdef EIGEN_VECTORIZE_AVX EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("AVX") != std::string::npos) - << "Not using AVX instructions:" << Eigen::SimdInstructionSetsInUse(); -#ifndef EIGEN_VECTORIZE_AVX - EXPECT_TRUE(false) - << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse(); +#else + EXPECT_TRUE(true) #endif + << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse(); +} + +TEST(flags, NEON) { +#ifdef EIGEN_VECTORIZE_NEON + EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("NEON") != std::string::npos) +#else + EXPECT_TRUE(true) +#endif + << "Not using EIGEN_VECTORIZE_NEON:" << Eigen::SimdInstructionSetsInUse(); } #define QUOTE(...) STRFY(__VA_ARGS__) diff --git a/src/include/cppduals/thirdparty/CMakeLists.txt b/src/include/cppduals/thirdparty/CMakeLists.txt index 12b274152..45cd32bc3 100644 --- a/src/include/cppduals/thirdparty/CMakeLists.txt +++ b/src/include/cppduals/thirdparty/CMakeLists.txt @@ -4,12 +4,19 @@ # 3.10 adds support for "gtest_discover_tests" which enumerates the tests inside # of the code and adds them to ctest. # -cmake_minimum_required (VERSION 3.10) +cmake_minimum_required (VERSION 3.14) project (cppduals_thirdparty) include (ExternalProject) get_directory_property (hasParent PARENT_DIRECTORY) +if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.23) + cmake_policy (SET CMP0135 NEW) +endif () +if (CMAKE_XCODE_BUILD_SYSTEM VERSION_GREATER_EQUAL 12) + cmake_policy (SET CMP0114 NEW) +endif () + set (DEPS_ROOT "${CMAKE_BINARY_DIR}/root") if (hasParent) set (DEPS_ROOT "${CMAKE_BINARY_DIR}/thirdparty/root" PARENT_SCOPE) @@ -21,98 +28,109 @@ else (NOT WIN32) set (DOWNLOAD_DIR "C:/Downloads") endif (NOT WIN32) +include(FetchContent) + # -# Google test (https://github.com/google/googletest/blob/master/googletest/README.md) -# +# Google test +# https://google.github.io/googletest/quickstart-cmake.html -# Download and unpack googletest at configure time -configure_file (CMakeLists-gt.txt.in googletest-download/CMakeLists.txt) -execute_process (COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" . - RESULT_VARIABLE result - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/googletest-download) -if (result) - message (FATAL_ERROR "CMake step for googletest failed: ${result}") -endif () -execute_process (COMMAND ${CMAKE_COMMAND} --build . - RESULT_VARIABLE result - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/googletest-download) -if (result) - message (FATAL_ERROR "Build step for googletest failed: ${result}") -endif () +# +# Configure google-test +FetchContent_Declare( + googletest + URL https://github.com/google/googletest/archive/refs/tags/v1.15.2.zip + URL_HASH MD5=eb1c5c237d13ed12bf492d3997ca6b0d + DOWNLOAD_NAME googletest-v1.15.2.zip + DOWNLOAD_DIR "$ENV{HOME}/Downloads" +) -# Prevent overriding the parent project's compiler/linker -# settings on Windows -set (gtest_force_shared_crt ON CACHE BOOL "" FORCE) +# For Windows: Prevent overriding the parent project's compiler/linker settings +set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) +FetchContent_MakeAvailable(googletest) -# Add googletest directly to our build. This defines -# the gtest and gtest_main targets. -add_subdirectory ( - ${CMAKE_CURRENT_BINARY_DIR}/googletest-src - ${CMAKE_CURRENT_BINARY_DIR}/googletest-build - EXCLUDE_FROM_ALL) - -# The gtest/gtest_main targets carry header search path -# dependencies automatically when using CMake 2.8.11 or -# later. Otherwise we have to add them here ourselves. -if (CMAKE_VERSION VERSION_LESS 2.8.11) - include_directories ("${gtest_SOURCE_DIR}/include") -endif () - -# Can simply link against gtest or gtest_main as needed. Eg -#add_executable (example example.cpp) -#target_link_libraries (example gtest_main) -#add_test (NAME example_test COMMAND example) +# Enable testing +include(GoogleTest) # # Eigen # -if (CPPDUALS_EIGEN_LATEST) - set (EIGEN_URL http://bitbucket.org/eigen/eigen/get/default.tar.bz2) - #set (EIGEN_MD5 ffc83130dcd37b694c6cf7e905099af9) +if (TRUE) + set (EIGEN_URL https://gitlab.com/libeigen/eigen/-/archive/3.3.8/eigen-3.3.8.tar.bz2) + set (EIGEN_MD5 432ef01499d514f4606343276afa0ec3) + set (EIGEN_MAX_CXX 17) else () - set (EIGEN_URL http://bitbucket.org/eigen/eigen/get/3.3.7.tar.bz2) - set (EIGEN_MD5 05b1f7511c93980c385ebe11bd3c93fa) + set (EIGEN_URL https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2) + set (EIGEN_MD5 b9e98a200d2455f06db9c661c5610496) + set (EIGEN_MAX_CXX 17) endif () - +# +# Eigen +# ExternalProject_Add (eigenX PREFIX eigenX URL ${EIGEN_URL} - #URL_HASH MD5=${EIGEN_MD5} + URL_HASH MD5=${EIGEN_MD5} DOWNLOAD_DIR "$ENV{HOME}/Downloads" CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + LOG_DOWNLOAD ON ) ExternalProject_Get_Property (eigenX source_dir) -if (hasParent AND NOT EIGEN3_INCLUDE_DIRS) - set (EIGEN3_INCLUDE_DIRS "${source_dir}" PARENT_SCOPE) +if (true) # || hasParent + add_library (eigen INTERFACE IMPORTED GLOBAL) + add_dependencies (eigen eigenX) + set_property (TARGET eigen APPEND PROPERTY INTERFACE_INCLUDE_DIRECTORIES "${source_dir}") + if (XCODE) + set (IOFORMAT "IOFormat\(FullPrecision, DontAlignCols, \", \", \"\\$\\\\n\", \"\", \"\", \"[\", \"]\"\)") + else () + set (IOFORMAT "IOFormat\(FullPrecision, DontAlignCols, \", \", \"\\$\\n\", \"\", \"\", \"[\", \"]\"\)") + endif () + #target_compile_definitions (eigen INTERFACE EIGEN_DEFAULT_IO_FORMAT=${IOFORMAT}) + #target_compile_definitions (eigen INTERFACE EIGEN_DEFAULT_IO_FORMAT=EIGEN_IO_FORMAT_OCTAVE) + set_property (TARGET eigen APPEND PROPERTY + INTERFACE_COMPILE_DEFINITIONS EIGEN_DEFAULT_IO_FORMAT=${IOFORMAT}) +endif () + +# if c++20, disable warning -Wdeprecated-enum-enum-conversion to eigen +if (CMAKE_CXX_STANDARD GREATER_EQUAL 20) + target_compile_options (eigen INTERFACE -Wno-deprecated-enum-enum-conversion) endif () # -# Eigen-Expokit +# expokit # -set (EEX_SHA 72bf6e445d5ae84218dcbd74580720491e0074db ) +#set (EEX_SHA ee28baa3bf29561501e17e5c68c2e54c85daae19 ) newer? used by spindropsSDL. md5=cebd15f9b5068c0e327753244ff6d394 +set (EEX_SHA c157dec0057be6e183a1ea2a5de353fac7e5e3a7 ) +set (EEX_MD5 89484e51f706398284235b96bc805515 ) ExternalProject_Add (expokitX PREFIX expokitX - URL https://gitlab.com/api/v4/projects/tesch1%2Feigen-expokit/repository/archive.tbz2?sha=${EEX_SHA} - #URL_HASH MD5=96b79de1d01547f6d658865b7caa02ee + URL https://gitlab.com/tesch1/eigen-expokit/-/archive/${EEX_SHA}/eigen-expokit.tar.bz2 + URL_HASH MD5=${EEX_MD5} + DOWNLOAD_NAME eigen-expokit-${EEX_SHA}.tar.bz2 DOWNLOAD_DIR "$ENV{HOME}/Downloads" + UPDATE_COMMAND "" + PATCH_COMMAND "" CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" + LOG_DOWNLOAD ON ) ExternalProject_Get_Property (expokitX source_dir) if (hasParent) - set (EXPOKIT_INCLUDE_DIR "${source_dir}" PARENT_SCOPE) -endif() + add_library (expokit INTERFACE IMPORTED GLOBAL) + add_dependencies (expokit expokitX) + set_property (TARGET expokit APPEND PROPERTY INTERFACE_INCLUDE_DIRECTORIES "${source_dir}") +endif (hasParent) + # # fmt # ExternalProject_Add (fmtX PREFIX fmtX - URL https://github.com/fmtlib/fmt/archive/6.1.1.tar.gz - URL_HASH MD5=acfb83d44afdca171ee26c597c931e7c + URL https://github.com/fmtlib/fmt/archive/11.1.4.tar.gz + URL_HASH MD5=10c2ae163accd3b82e6b8b4dff877645 DOWNLOAD_DIR ${DOWNLOAD_DIR} CONFIGURE_COMMAND "" BUILD_COMMAND "" @@ -121,7 +139,6 @@ ExternalProject_Add (fmtX ExternalProject_Get_Property (fmtX source_dir) ExternalProject_Get_Property (fmtX binary_dir) if (hasParent) - message (" FMT3_INCLUDE_DIRS: ${source_dir}") add_subdirectory (${source_dir} ${binary_dir} EXCLUDE_FROM_ALL) endif () @@ -131,19 +148,20 @@ if (CPPDUALS_BENCHMARK) # ExternalProject_Add (benchmarkX PREFIX benchmarkX - URL "http://github.com/google/benchmark/archive/v1.5.0.tar.gz" - URL_HASH MD5=eb1466370f3ae31e74557baa29729e9e + URL "http://github.com/google/benchmark/archive/v1.9.1.tar.gz" + URL_HASH MD5=92000ef8b4a7b1e9229972f8943070a7 DOWNLOAD_DIR ${DOWNLOAD_DIR} - CMAKE_ARGS --target install -DBENCHMARK_ENABLE_GTEST_TESTS=OFF -DCMAKE_BUILD_TYPE=Release -DBENCHMARK_USE_LIBCXX=${CPPDUALS_USE_LIBCXX} - "-DCMAKE_INSTALL_PREFIX=" - INSTALL_DIR "${DEPS_ROOT}" + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" ) ExternalProject_Get_Property (benchmarkX source_dir) - ExternalProject_Get_Property (benchmarkX install_dir) + ExternalProject_Get_Property (benchmarkX binary_dir) if (hasParent) - set (BENCHMARK_SRC_DIR "${source_dir}" PARENT_SCOPE) - set (BENCHMARK_INC_DIR "${install_dir}/include" PARENT_SCOPE) - message (" BENCHMARK_SRC_DIR: ${BENCHMARK_SRC_DIR}") + # https://github.com/google/benchmark#requirements + set (BENCHMARK_ENABLE_GTEST_TESTS OFF) + set (BENCHMARK_USE_LIBCXX ${CPPDUALS_USE_LIBCXX}) + add_subdirectory (${source_dir} ${binary_dir} EXCLUDE_FROM_ALL) endif() if (Boost_FOUND AND NO) @@ -158,18 +176,6 @@ if (CPPDUALS_BENCHMARK) set (Boost_INCLUDE_DIRS ${Boost_INCLUDE_DIRS} PARENT_SCOPE) endif () - # piranha - ExternalProject_Add (piranhaX PREFIX piranhaX - URL https://github.com/bluescarni/piranha/archive/v0.11.tar.gz - URL_HASH MD5=33482f719f6b8a6a5316f9bd148f5b10 - DOWNLOAD_DIR "$ENV{HOME}/Downloads" - CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" - ) - ExternalProject_Get_Property (piranhaX source_dir) - if (hasParent) - set (PIRANHA_INCLUDE_DIR "${source_dir}/include" PARENT_SCOPE) - endif () - # AuDi ExternalProject_Add (audiX PREFIX audiX URL https://github.com/darioizzo/audi/archive/v1.6.5.tar.gz