diff --git a/cmake/configure.cmake b/cmake/configure.cmake index 1814656d244088bbf06611109d18caea73ad6a45..8200d48cbb3947e723abe09e511b73312191887a 100644 --- a/cmake/configure.cmake +++ b/cmake/configure.cmake @@ -70,10 +70,6 @@ endif() if(WITH_GPU) add_definitions(-DPADDLE_WITH_CUDA) add_definitions(-DEIGEN_USE_GPU) - # The compiler fully support const expressions since c++14, - # but Eigen use some const expressions such as std::max and std::min, which are not supported in c++11 - # use following definition to set EIGEN_HAS_CONSTEXPR=0 to avoid compilation error in c++11 - add_definitions(-DEIGEN_MAX_CPP_VER=11) FIND_PACKAGE(CUDA REQUIRED) diff --git a/cmake/external/eigen.cmake b/cmake/external/eigen.cmake index 30e74a3c241dcb48aa3bb9356dcc42a46283045c..631803da31d5a0fad1534d266f70f4708882f7e8 100644 --- a/cmake/external/eigen.cmake +++ b/cmake/external/eigen.cmake @@ -49,9 +49,14 @@ elseif(LINUX) # refer to: https://gitlab.com/libeigen/eigen/-/blob/4da2c6b1974827b1999bab652a3d4703e1992d26/Eigen/src/Core/arch/SSE/PacketMath.h#L33-60 # add -fabi-version=4 could avoid above error, but will cause "double free corruption" when compile with gcc8 # so use following patch to solve compilation error with different version of gcc. - file(TO_NATIVE_PATH ${PADDLE_SOURCE_DIR}/patches/eigen/Geometry_SSE.h native_src) - file(TO_NATIVE_PATH ${EIGEN_SOURCE_DIR}/Eigen/src/Geometry/arch/Geometry_SSE.h native_dst) - set(EIGEN_PATCH_COMMAND cp ${native_src} ${native_dst}) + file(TO_NATIVE_PATH ${PADDLE_SOURCE_DIR}/patches/eigen/Geometry_SSE.h native_src1) + file(TO_NATIVE_PATH ${EIGEN_SOURCE_DIR}/Eigen/src/Geometry/arch/Geometry_SSE.h native_dst1) + # The compiler fully support const expressions since c++14, + # but Eigen use some const expressions such as std::max and std::min, which are not supported in c++11 + # add patch to avoid compilation error in c++11 + file(TO_NATIVE_PATH ${PADDLE_SOURCE_DIR}/patches/eigen/MathFunctions.h native_src2) + file(TO_NATIVE_PATH ${EIGEN_SOURCE_DIR}/Eigen/src/Core/MathFunctions.h native_dst2) + set(EIGEN_PATCH_COMMAND cp ${native_src1} ${native_dst1} && cp ${native_src2} ${native_dst2}) endif() set(EIGEN_INCLUDE_DIR ${EIGEN_SOURCE_DIR}) diff --git a/patches/eigen/MathFunctions.h b/patches/eigen/MathFunctions.h new file mode 100644 index 0000000000000000000000000000000000000000..9f6a4d0e3328ffbe864273b69dc0033d50b813ae --- /dev/null +++ b/patches/eigen/MathFunctions.h @@ -0,0 +1,1938 @@ +// Copyright (c) 2020 PaddlePaddle Authors. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2006-2010 Benoit Jacob +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_MATHFUNCTIONS_H +#define EIGEN_MATHFUNCTIONS_H + +// source: http://www.geom.uiuc.edu/~huberty/math5337/groupe/digits.html +// TODO this should better be moved to NumTraits +#define EIGEN_PI \ + 3.141592653589793238462643383279502884197169399375105820974944592307816406L + +namespace Eigen { + +// On WINCE, std::abs is defined for int only, so let's defined our own +// overloads: +// This issue has been confirmed with MSVC 2008 only, but the issue might exist +// for more recent versions too. +#if EIGEN_OS_WINCE && EIGEN_COMP_MSVC && EIGEN_COMP_MSVC <= 1500 +long abs(long x) { return (labs(x)); } +double abs(double x) { return (fabs(x)); } +float abs(float x) { return (fabsf(x)); } +long double abs(long double x) { return (fabsl(x)); } +#endif + +namespace internal { + +/** \internal \class global_math_functions_filtering_base + * + * What it does: + * Defines a typedef 'type' as follows: + * - if type T has a member typedef + * Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl, then + * global_math_functions_filtering_base::type is a typedef for it. + * - otherwise, global_math_functions_filtering_base::type is a typedef for + * T. + * + * How it's used: + * To allow to defined the global math functions (like sin...) in certain + * cases, like the Array expressions. + * When you do sin(array1+array2), the object array1+array2 has a complicated + * expression type, all what you want to know + * is that it inherits ArrayBase. So we implement a partial specialization of + * sin_impl for ArrayBase. + * So we must make sure to use sin_impl > and not + * sin_impl, otherwise our partial specialization + * won't be used. How does sin know that? That's exactly what + * global_math_functions_filtering_base tells it. + * + * How it's implemented: + * SFINAE in the style of enable_if. Highly susceptible of breaking compilers. + * With GCC, it sure does work, but if you replace + * the typename dummy by an integer template parameter, it doesn't work + * anymore! + */ + +template +struct global_math_functions_filtering_base { + typedef T type; +}; + +template +struct always_void { + typedef void type; +}; + +template +struct global_math_functions_filtering_base< + T, + typename always_void< + typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl>:: + type> { + typedef typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl type; +}; + +#define EIGEN_MATHFUNC_IMPL(func, scalar) \ + Eigen::internal::func##_impl< \ + typename Eigen::internal::global_math_functions_filtering_base< \ + scalar>::type> +#define EIGEN_MATHFUNC_RETVAL(func, scalar) \ + typename Eigen::internal::func##_retval< \ + typename Eigen::internal::global_math_functions_filtering_base< \ + scalar>::type>::type + +/**************************************************************************** +* Implementation of real * +****************************************************************************/ + +template ::IsComplex> +struct real_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { return x; } +}; + +template +struct real_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + using std::real; + return real(x); + } +}; + +template +struct real_impl : real_default_impl {}; + +#if defined(EIGEN_GPU_COMPILE_PHASE) +template +struct real_impl> { + typedef T RealScalar; + EIGEN_DEVICE_FUNC + static inline T run(const std::complex& x) { return x.real(); } +}; +#endif + +template +struct real_retval { + typedef typename NumTraits::Real type; +}; + +/**************************************************************************** +* Implementation of imag * +****************************************************************************/ + +template ::IsComplex> +struct imag_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar&) { return RealScalar(0); } +}; + +template +struct imag_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + using std::imag; + return imag(x); + } +}; + +template +struct imag_impl : imag_default_impl {}; + +#if defined(EIGEN_GPU_COMPILE_PHASE) +template +struct imag_impl> { + typedef T RealScalar; + EIGEN_DEVICE_FUNC + static inline T run(const std::complex& x) { return x.imag(); } +}; +#endif + +template +struct imag_retval { + typedef typename NumTraits::Real type; +}; + +/**************************************************************************** +* Implementation of real_ref * +****************************************************************************/ + +template +struct real_ref_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar& run(Scalar& x) { + return reinterpret_cast(&x)[0]; + } + EIGEN_DEVICE_FUNC + static inline const RealScalar& run(const Scalar& x) { + return reinterpret_cast(&x)[0]; + } +}; + +template +struct real_ref_retval { + typedef typename NumTraits::Real& type; +}; + +/**************************************************************************** +* Implementation of imag_ref * +****************************************************************************/ + +template +struct imag_ref_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar& run(Scalar& x) { + return reinterpret_cast(&x)[1]; + } + EIGEN_DEVICE_FUNC + static inline const RealScalar& run(const Scalar& x) { + return reinterpret_cast(&x)[1]; + } +}; + +template +struct imag_ref_default_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(Scalar&) { return Scalar(0); } + EIGEN_DEVICE_FUNC + static inline const Scalar run(const Scalar&) { return Scalar(0); } +}; + +template +struct imag_ref_impl + : imag_ref_default_impl::IsComplex> {}; + +template +struct imag_ref_retval { + typedef typename NumTraits::Real& type; +}; + +/**************************************************************************** +* Implementation of conj * +****************************************************************************/ + +template ::IsComplex> +struct conj_default_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { return x; } +}; + +template +struct conj_default_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { + using std::conj; + return conj(x); + } +}; + +template +struct conj_impl : conj_default_impl {}; + +#if defined(EIGEN_GPU_COMPILE_PHASE) +template +struct conj_impl> { + EIGEN_DEVICE_FUNC + static inline std::complex run(const std::complex& x) { + return std::complex(x.real(), -x.imag()); + } +}; +#endif + +template +struct conj_retval { + typedef Scalar type; +}; + +/**************************************************************************** +* Implementation of abs2 * +****************************************************************************/ + +template +struct abs2_impl_default { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { return x * x; } +}; + +template +struct abs2_impl_default // IsComplex +{ + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + return x.real() * x.real() + x.imag() * x.imag(); + } +}; + +template +struct abs2_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + return abs2_impl_default::IsComplex>::run(x); + } +}; + +template +struct abs2_retval { + typedef typename NumTraits::Real type; +}; + +/**************************************************************************** +* Implementation of norm1 * +****************************************************************************/ + +template +struct norm1_default_impl; + +template +struct norm1_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + EIGEN_USING_STD_MATH(abs); + return abs(x.real()) + abs(x.imag()); + } +}; + +template +struct norm1_default_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { + EIGEN_USING_STD_MATH(abs); + return abs(x); + } +}; + +template +struct norm1_impl : norm1_default_impl::IsComplex> {}; + +template +struct norm1_retval { + typedef typename NumTraits::Real type; +}; + +/**************************************************************************** +* Implementation of hypot * +****************************************************************************/ + +template +struct hypot_impl; + +template +struct hypot_retval { + typedef typename NumTraits::Real type; +}; + +/**************************************************************************** +* Implementation of cast * +****************************************************************************/ + +template +struct cast_impl { + EIGEN_DEVICE_FUNC + static inline NewType run(const OldType& x) { + return static_cast(x); + } +}; + +// here, for once, we're plainly returning NewType: we don't want cast to do +// weird things. + +template +EIGEN_DEVICE_FUNC inline NewType cast(const OldType& x) { + return cast_impl::run(x); +} + +/**************************************************************************** +* Implementation of round * +****************************************************************************/ + +#if EIGEN_HAS_CXX11_MATH +template +struct round_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { + EIGEN_STATIC_ASSERT((!NumTraits::IsComplex), + NUMERIC_TYPE_MUST_BE_REAL) + EIGEN_USING_STD_MATH(round); + return round(x); + } +}; +#else +template +struct round_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { + EIGEN_STATIC_ASSERT((!NumTraits::IsComplex), + NUMERIC_TYPE_MUST_BE_REAL) + EIGEN_USING_STD_MATH(floor); + EIGEN_USING_STD_MATH(ceil); + return (x > Scalar(0)) ? floor(x + Scalar(0.5)) : ceil(x - Scalar(0.5)); + } +}; +#endif + +template +struct round_retval { + typedef Scalar type; +}; + +/**************************************************************************** +* Implementation of rint * +****************************************************************************/ + +template +struct rint_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { + EIGEN_STATIC_ASSERT((!NumTraits::IsComplex), + NUMERIC_TYPE_MUST_BE_REAL) +#if EIGEN_HAS_CXX11_MATH + EIGEN_USING_STD_MATH(rint); +#endif + return rint(x); + } +}; + +#if !EIGEN_HAS_CXX11_MATH +template <> +struct rint_impl { + EIGEN_DEVICE_FUNC + static inline double run(const double& x) { return ::rint(x); } +}; +template <> +struct rint_impl { + EIGEN_DEVICE_FUNC + static inline float run(const float& x) { return ::rintf(x); } +}; +#endif + +template +struct rint_retval { + typedef Scalar type; +}; + +/**************************************************************************** +* Implementation of arg * +****************************************************************************/ + +#if EIGEN_HAS_CXX11_MATH +template +struct arg_impl { + EIGEN_DEVICE_FUNC + static inline Scalar run(const Scalar& x) { +#if defined(EIGEN_HIP_DEVICE_COMPILE) + // HIP does not seem to have a native device side implementation for the + // math routine "arg" + using std::arg; +#else + EIGEN_USING_STD_MATH(arg); +#endif + return arg(x); + } +}; +#else +template ::IsComplex> +struct arg_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + return (x < Scalar(0)) ? Scalar(EIGEN_PI) : Scalar(0); + } +}; + +template +struct arg_default_impl { + typedef typename NumTraits::Real RealScalar; + EIGEN_DEVICE_FUNC + static inline RealScalar run(const Scalar& x) { + EIGEN_USING_STD_MATH(arg); + return arg(x); + } +}; + +template +struct arg_impl : arg_default_impl {}; +#endif + +template +struct arg_retval { + typedef typename NumTraits::Real type; +}; + +/**************************************************************************** +* Implementation of expm1 * +****************************************************************************/ + +// This implementation is based on GSL Math's expm1. +namespace std_fallback { +// fallback expm1 implementation in case there is no expm1(Scalar) function in +// namespace of Scalar, +// or that there is no suitable std::expm1 function available. Implementation +// attributed to Kahan. See: http://www.plunk.org/~hatch/rightway.php. +template +EIGEN_DEVICE_FUNC inline Scalar expm1(const Scalar& x) { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) + typedef typename NumTraits::Real RealScalar; + + EIGEN_USING_STD_MATH(exp); + Scalar u = exp(x); + if (numext::equal_strict(u, Scalar(1))) { + return x; + } + Scalar um1 = u - RealScalar(1); + if (numext::equal_strict(um1, Scalar(-1))) { + return RealScalar(-1); + } + + EIGEN_USING_STD_MATH(log); + Scalar logu = log(u); + return numext::equal_strict(u, logu) ? u : (u - RealScalar(1)) * x / logu; +} +} + +template +struct expm1_impl { + EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& x) { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) +#if EIGEN_HAS_CXX11_MATH + using std::expm1; +#else + using std_fallback::expm1; +#endif + return expm1(x); + } +}; + +// Specialization for complex types that are not supported by std::expm1. +template +struct expm1_impl> { + EIGEN_DEVICE_FUNC static inline std::complex run( + const std::complex& x) { + EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar) + RealScalar xr = x.real(); + RealScalar xi = x.imag(); + // expm1(z) = exp(z) - 1 + // = exp(x + i * y) - 1 + // = exp(x) * (cos(y) + i * sin(y)) - 1 + // = exp(x) * cos(y) - 1 + i * exp(x) * sin(y) + // Imag(expm1(z)) = exp(x) * sin(y) + // Real(expm1(z)) = exp(x) * cos(y) - 1 + // = exp(x) * cos(y) - 1. + // = expm1(x) + exp(x) * (cos(y) - 1) + // = expm1(x) + exp(x) * (2 * sin(y / 2) ** 2) + + // TODO better use numext::expm1 and numext::sin (but that would require + // forward declarations or moving this specialization down). + RealScalar erm1 = expm1_impl::run(xr); + RealScalar er = erm1 + RealScalar(1.); + EIGEN_USING_STD_MATH(sin); + RealScalar sin2 = sin(xi / RealScalar(2.)); + sin2 = sin2 * sin2; + RealScalar s = sin(xi); + RealScalar real_part = erm1 - RealScalar(2.) * er * sin2; + return std::complex(real_part, er * s); + } +}; + +template +struct expm1_retval { + typedef Scalar type; +}; + +/**************************************************************************** +* Implementation of log1p * +****************************************************************************/ + +namespace std_fallback { +// fallback log1p implementation in case there is no log1p(Scalar) function in +// namespace of Scalar, +// or that there is no suitable std::log1p function available +template +EIGEN_DEVICE_FUNC inline Scalar log1p(const Scalar& x) { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) + typedef typename NumTraits::Real RealScalar; + EIGEN_USING_STD_MATH(log); + Scalar x1p = RealScalar(1) + x; + Scalar log_1p = log(x1p); + const bool is_small = numext::equal_strict(x1p, Scalar(1)); + const bool is_inf = numext::equal_strict(x1p, log_1p); + return (is_small || is_inf) ? x : x * (log_1p / (x1p - RealScalar(1))); +} +} + +template +struct log1p_impl { + EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& x) { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) +#if EIGEN_HAS_CXX11_MATH + using std::log1p; +#else + using std_fallback::log1p; +#endif + return log1p(x); + } +}; + +// Specialization for complex types that are not supported by std::log1p. +template +struct log1p_impl> { + EIGEN_DEVICE_FUNC static inline std::complex run( + const std::complex& x) { + EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar) + return std_fallback::log1p(x); + } +}; + +template +struct log1p_retval { + typedef Scalar type; +}; + +/**************************************************************************** +* Implementation of pow * +****************************************************************************/ + +template ::IsInteger&& NumTraits::IsInteger> +struct pow_impl { + // typedef Scalar retval; + typedef typename ScalarBinaryOpTraits< + ScalarX, + ScalarY, + internal::scalar_pow_op>::ReturnType result_type; + static EIGEN_DEVICE_FUNC inline result_type run(const ScalarX& x, + const ScalarY& y) { + EIGEN_USING_STD_MATH(pow); + return pow(x, y); + } +}; + +template +struct pow_impl { + typedef ScalarX result_type; + static EIGEN_DEVICE_FUNC inline ScalarX run(ScalarX x, ScalarY y) { + ScalarX res(1); + eigen_assert(!NumTraits::IsSigned || y >= 0); + if (y & 1) res *= x; + y >>= 1; + while (y) { + x *= x; + if (y & 1) res *= x; + y >>= 1; + } + return res; + } +}; + +/**************************************************************************** +* Implementation of random * +****************************************************************************/ + +template +struct random_default_impl {}; + +template +struct random_impl : random_default_impl::IsComplex, + NumTraits::IsInteger> {}; + +template +struct random_retval { + typedef Scalar type; +}; + +template +inline EIGEN_MATHFUNC_RETVAL(random, Scalar) + random(const Scalar& x, const Scalar& y); +template +inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(); + +template +struct random_default_impl { + static inline Scalar run(const Scalar& x, const Scalar& y) { + return x + (y - x) * Scalar(std::rand()) / Scalar(RAND_MAX); + } + static inline Scalar run() { + return run(Scalar(NumTraits::IsSigned ? -1 : 0), Scalar(1)); + } +}; + +enum { + meta_floor_log2_terminate, + meta_floor_log2_move_up, + meta_floor_log2_move_down, + meta_floor_log2_bogus +}; + +template +struct meta_floor_log2_selector { + enum { + middle = (lower + upper) / 2, + value = (upper <= lower + 1) + ? int(meta_floor_log2_terminate) + : (n < (1 << middle)) ? int(meta_floor_log2_move_down) + : (n == 0) ? int(meta_floor_log2_bogus) + : int(meta_floor_log2_move_up) + }; +}; + +template ::value> +struct meta_floor_log2 {}; + +template +struct meta_floor_log2 { + enum { + value = meta_floor_log2< + n, + lower, + meta_floor_log2_selector::middle>::value + }; +}; + +template +struct meta_floor_log2 { + enum { + value = meta_floor_log2::middle, + upper>::value + }; +}; + +template +struct meta_floor_log2 { + enum { + value = (n >= ((unsigned int)(1) << (lower + 1))) ? lower + 1 : lower + }; +}; + +template +struct meta_floor_log2 { + // no value, error at compile time +}; + +template +struct random_default_impl { + static inline Scalar run(const Scalar& x, const Scalar& y) { + if (y <= x) return x; + // ScalarU is the unsigned counterpart of Scalar, possibly Scalar itself. + typedef typename make_unsigned::type ScalarU; + // ScalarX is the widest of ScalarU and unsigned int. + // We'll deal only with ScalarX and unsigned int below thus avoiding signed + // types and arithmetic and signed overflows (which are undefined behavior). + typedef typename conditional<(ScalarU(-1) > unsigned(-1)), + ScalarU, + unsigned>::type ScalarX; + // The following difference doesn't overflow, provided our integer types are + // two's + // complement and have the same number of padding bits in signed and + // unsigned variants. + // This is the case in most modern implementations of C++. + ScalarX range = ScalarX(y) - ScalarX(x); + ScalarX offset = 0; + ScalarX divisor = 1; + ScalarX multiplier = 1; + const unsigned rand_max = RAND_MAX; + if (range <= rand_max) + divisor = (rand_max + 1) / (range + 1); + else + multiplier = 1 + range / (rand_max + 1); + // Rejection sampling. + do { + offset = (unsigned(std::rand()) * multiplier) / divisor; + } while (offset > range); + return Scalar(ScalarX(x) + offset); + } + + static inline Scalar run() { +#ifdef EIGEN_MAKING_DOCS + return run(Scalar(NumTraits::IsSigned ? -10 : 0), Scalar(10)); +#else + enum { + rand_bits = meta_floor_log2<(unsigned int)(RAND_MAX) + 1>::value, + scalar_bits = sizeof(Scalar) * CHAR_BIT, + shift = EIGEN_PLAIN_ENUM_MAX(0, int(rand_bits) - int(scalar_bits)), + offset = NumTraits::IsSigned + ? (1 << (EIGEN_PLAIN_ENUM_MIN(rand_bits, scalar_bits) - 1)) + : 0}; + return Scalar((std::rand() >> shift) - offset); +#endif + } +}; + +template +struct random_default_impl { + static inline Scalar run(const Scalar& x, const Scalar& y) { + return Scalar(random(x.real(), y.real()), random(x.imag(), y.imag())); + } + static inline Scalar run() { + typedef typename NumTraits::Real RealScalar; + return Scalar(random(), random()); + } +}; + +template +inline EIGEN_MATHFUNC_RETVAL(random, Scalar) + random(const Scalar& x, const Scalar& y) { + return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y); +} + +template +inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random() { + return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(); +} + +// Implementation of is* functions + +// std::is* do not work with fast-math and gcc, std::is* are available on MSVC +// 2013 and newer, as well as in clang. +#if (EIGEN_HAS_CXX11_MATH && \ + !(EIGEN_COMP_GNUC_STRICT && __FINITE_MATH_ONLY__)) || \ + (EIGEN_COMP_MSVC >= 1800) || (EIGEN_COMP_CLANG) +#define EIGEN_USE_STD_FPCLASSIFY 1 +#else +#define EIGEN_USE_STD_FPCLASSIFY 0 +#endif + +template +EIGEN_DEVICE_FUNC + typename internal::enable_if::value, bool>::type + isnan_impl(const T&) { + return false; +} + +template +EIGEN_DEVICE_FUNC + typename internal::enable_if::value, bool>::type + isinf_impl(const T&) { + return false; +} + +template +EIGEN_DEVICE_FUNC + typename internal::enable_if::value, bool>::type + isfinite_impl(const T&) { + return true; +} + +template +EIGEN_DEVICE_FUNC + typename internal::enable_if<(!internal::is_integral::value) && + (!NumTraits::IsComplex), + bool>::type + isfinite_impl(const T& x) { +#if defined(EIGEN_GPU_COMPILE_PHASE) + return (::isfinite)(x); +#elif EIGEN_USE_STD_FPCLASSIFY + using std::isfinite; + return isfinite EIGEN_NOT_A_MACRO(x); +#else + return x <= NumTraits::highest() && x >= NumTraits::lowest(); +#endif +} + +template +EIGEN_DEVICE_FUNC + typename internal::enable_if<(!internal::is_integral::value) && + (!NumTraits::IsComplex), + bool>::type + isinf_impl(const T& x) { +#if defined(EIGEN_GPU_COMPILE_PHASE) + return (::isinf)(x); +#elif EIGEN_USE_STD_FPCLASSIFY + using std::isinf; + return isinf EIGEN_NOT_A_MACRO(x); +#else + return x > NumTraits::highest() || x < NumTraits::lowest(); +#endif +} + +template +EIGEN_DEVICE_FUNC + typename internal::enable_if<(!internal::is_integral::value) && + (!NumTraits::IsComplex), + bool>::type + isnan_impl(const T& x) { +#if defined(EIGEN_GPU_COMPILE_PHASE) + return (::isnan)(x); +#elif EIGEN_USE_STD_FPCLASSIFY + using std::isnan; + return isnan EIGEN_NOT_A_MACRO(x); +#else + return x != x; +#endif +} + +#if (!EIGEN_USE_STD_FPCLASSIFY) + +#if EIGEN_COMP_MSVC + +template +EIGEN_DEVICE_FUNC bool isinf_msvc_helper(T x) { + return _fpclass(x) == _FPCLASS_NINF || _fpclass(x) == _FPCLASS_PINF; +} + +// MSVC defines a _isnan builtin function, but for double only +EIGEN_DEVICE_FUNC inline bool isnan_impl(const long double& x) { + return _isnan(x) != 0; +} +EIGEN_DEVICE_FUNC inline bool isnan_impl(const double& x) { + return _isnan(x) != 0; +} +EIGEN_DEVICE_FUNC inline bool isnan_impl(const float& x) { + return _isnan(x) != 0; +} + +EIGEN_DEVICE_FUNC inline bool isinf_impl(const long double& x) { + return isinf_msvc_helper(x); +} +EIGEN_DEVICE_FUNC inline bool isinf_impl(const double& x) { + return isinf_msvc_helper(x); +} +EIGEN_DEVICE_FUNC inline bool isinf_impl(const float& x) { + return isinf_msvc_helper(x); +} + +#elif (defined __FINITE_MATH_ONLY__ && __FINITE_MATH_ONLY__ && EIGEN_COMP_GNUC) + +#if EIGEN_GNUC_AT_LEAST(5, 0) +#define EIGEN_TMP_NOOPT_ATTRIB \ + EIGEN_DEVICE_FUNC inline __attribute__((optimize("no-finite-math-only"))) +#else +// NOTE the inline qualifier and noinline attribute are both needed: the former +// is to avoid linking issue (duplicate symbol), +// while the second prevent too aggressive optimizations in fast-math mode: +#define EIGEN_TMP_NOOPT_ATTRIB \ + EIGEN_DEVICE_FUNC inline \ + __attribute__((noinline, optimize("no-finite-math-only"))) +#endif + +template <> +EIGEN_TMP_NOOPT_ATTRIB bool isnan_impl(const long double& x) { + return __builtin_isnan(x); +} +template <> +EIGEN_TMP_NOOPT_ATTRIB bool isnan_impl(const double& x) { + return __builtin_isnan(x); +} +template <> +EIGEN_TMP_NOOPT_ATTRIB bool isnan_impl(const float& x) { + return __builtin_isnan(x); +} +template <> +EIGEN_TMP_NOOPT_ATTRIB bool isinf_impl(const double& x) { + return __builtin_isinf(x); +} +template <> +EIGEN_TMP_NOOPT_ATTRIB bool isinf_impl(const float& x) { + return __builtin_isinf(x); +} +template <> +EIGEN_TMP_NOOPT_ATTRIB bool isinf_impl(const long double& x) { + return __builtin_isinf(x); +} + +#undef EIGEN_TMP_NOOPT_ATTRIB + +#endif + +#endif + +// The following overload are defined at the end of this file +template +EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex& x); +template +EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex& x); +template +EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex& x); + +template +T generic_fast_tanh_float(const T& a_x); +} // end namespace internal + +/**************************************************************************** +* Generic math functions * +****************************************************************************/ + +namespace numext { + +#if (!defined(EIGEN_GPUCC)) +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T& x, const T& y) { + EIGEN_USING_STD_MATH(min); + return min EIGEN_NOT_A_MACRO(x, y); +} + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T& x, const T& y) { + EIGEN_USING_STD_MATH(max); + return max EIGEN_NOT_A_MACRO(x, y); +} +#else +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T& x, const T& y) { + return y < x ? y : x; +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float mini(const float& x, + const float& y) { + return fminf(x, y); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double mini(const double& x, + const double& y) { + return fmin(x, y); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE long double mini(const long double& x, + const long double& y) { +#if defined(EIGEN_HIPCC) + // no "fminl" on HIP yet + return (x < y) ? x : y; +#else + return fminl(x, y); +#endif +} + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T& x, const T& y) { + return x < y ? y : x; +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float maxi(const float& x, + const float& y) { + return fmaxf(x, y); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double maxi(const double& x, + const double& y) { + return fmax(x, y); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE long double maxi(const long double& x, + const long double& y) { +#if defined(EIGEN_HIPCC) + // no "fmaxl" on HIP yet + return (x > y) ? x : y; +#else + return fmaxl(x, y); +#endif +} +#endif + +#if defined(SYCL_DEVICE_ONLY) + +#define SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_char) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_short) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_int) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_long) +#define SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_char) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_short) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_int) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_long) +#define SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_uchar) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_ushort) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_uint) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_ulong) +#define SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_uchar) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_ushort) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_uint) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_ulong) +#define SYCL_SPECIALIZE_INTEGER_TYPES_BINARY(NAME, FUNC) \ + SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) \ + SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) +#define SYCL_SPECIALIZE_INTEGER_TYPES_UNARY(NAME, FUNC) \ + SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) \ + SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) +#define SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(NAME, FUNC) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_float) \ + SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_double) +#define SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(NAME, FUNC) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_float) \ + SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_double) +#define SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE( \ + NAME, FUNC, RET_TYPE) \ + SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, RET_TYPE, cl::sycl::cl_float) \ + SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, RET_TYPE, cl::sycl::cl_double) + +#define SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE) \ + template <> \ + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE RET_TYPE NAME(const ARG_TYPE& x) { \ + return cl::sycl::FUNC(x); \ + } + +#define SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, TYPE) \ + SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, TYPE, TYPE) + +#define SYCL_SPECIALIZE_GEN1_BINARY_FUNC( \ + NAME, FUNC, RET_TYPE, ARG_TYPE1, ARG_TYPE2) \ + template <> \ + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE RET_TYPE NAME(const ARG_TYPE1& x, \ + const ARG_TYPE2& y) { \ + return cl::sycl::FUNC(x, y); \ + } + +#define SYCL_SPECIALIZE_GEN2_BINARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE) \ + SYCL_SPECIALIZE_GEN1_BINARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE, ARG_TYPE) + +#define SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, TYPE) \ + SYCL_SPECIALIZE_GEN2_BINARY_FUNC(NAME, FUNC, TYPE, TYPE) + +SYCL_SPECIALIZE_INTEGER_TYPES_BINARY(mini, min) +SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(mini, fmin) +SYCL_SPECIALIZE_INTEGER_TYPES_BINARY(maxi, max) +SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(maxi, fmax) + +#endif + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(real, Scalar) + real(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline typename internal::add_const_on_value_type< + EIGEN_MATHFUNC_RETVAL(real_ref, Scalar)>::type +real_ref(const Scalar& x) { + return internal::real_ref_impl::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) + real_ref(Scalar& x) { + return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) + imag(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(arg, Scalar) + arg(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(arg, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline typename internal::add_const_on_value_type< + EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar)>::type +imag_ref(const Scalar& x) { + return internal::imag_ref_impl::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) + imag_ref(Scalar& x) { + return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) + conj(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) + abs2(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x); +} + +EIGEN_DEVICE_FUNC +inline bool abs2(bool x) { return x; } + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T absdiff(const T& x, const T& y) { + return x > y ? x - y : y - x; +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float absdiff(const float& x, + const float& y) { + return fabsf(x - y); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double absdiff(const double& x, + const double& y) { + return fabs(x - y); +} +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE long double absdiff( + const long double& x, const long double& y) { +#if defined(EIGEN_HIPCC) + // no "fabsl" on HIP yet + return (x > y) ? x : y; +#else + return fabsl(x - y); +#endif +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) + norm1(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) + hypot(const Scalar& x, const Scalar& y) { + return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(hypot, hypot) +#endif + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(log1p, Scalar) + log1p(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(log1p, Scalar)::run(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(log1p, log1p) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float log1p(const float& x) { + return ::log1pf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double log1p(const double& x) { + return ::log1p(x); +} +#endif + +template +EIGEN_DEVICE_FUNC inline + typename internal::pow_impl::result_type + pow(const ScalarX& x, const ScalarY& y) { + return internal::pow_impl::run(x, y); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(pow, pow) +#endif + +template +EIGEN_DEVICE_FUNC bool(isnan)(const T& x) { + return internal::isnan_impl(x); +} +template +EIGEN_DEVICE_FUNC bool(isinf)(const T& x) { + return internal::isinf_impl(x); +} +template +EIGEN_DEVICE_FUNC bool(isfinite)(const T& x) { + return internal::isfinite_impl(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(isnan, isnan, bool) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(isinf, isinf, bool) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(isfinite, isfinite, bool) +#endif + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(rint, Scalar) + rint(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(rint, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(round, Scalar) + round(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(round, Scalar)::run(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(round, round) +#endif + +template +EIGEN_DEVICE_FUNC T(floor)(const T& x) { + EIGEN_USING_STD_MATH(floor); + return floor(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(floor, floor) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float floor(const float& x) { + return ::floorf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double floor(const double& x) { + return ::floor(x); +} +#endif + +template +EIGEN_DEVICE_FUNC T(ceil)(const T& x) { + EIGEN_USING_STD_MATH(ceil); + return ceil(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(ceil, ceil) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float ceil(const float& x) { + return ::ceilf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double ceil(const double& x) { + return ::ceil(x); +} +#endif + +/** Log base 2 for 32 bits positive integers. + * Conveniently returns 0 for x==0. */ +inline int log2(int x) { + eigen_assert(x >= 0); + unsigned int v(x); + static const int table[32] = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, + 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, + 24, 7, 19, 27, 23, 6, 26, 5, 4, 31}; + v |= v >> 1; + v |= v >> 2; + v |= v >> 4; + v |= v >> 8; + v |= v >> 16; + return table[(v * 0x07C4ACDDU) >> 27]; +} + +/** \returns the square root of \a x. + * + * It is essentially equivalent to + * \code using std::sqrt; return sqrt(x); \endcode + * but slightly faster for float/double and some compilers (e.g., gcc), thanks + * to + * specializations when SSE is enabled. + * + * It's usage is justified in performance critical functions, like + * norm/normalize. + */ +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T sqrt(const T& x) { + EIGEN_USING_STD_MATH(sqrt); + return sqrt(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(sqrt, sqrt) +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T log(const T& x) { + EIGEN_USING_STD_MATH(log); + return log(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(log, log) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float log(const float& x) { + return ::logf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double log(const double& x) { + return ::log(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE + typename internal::enable_if::IsSigned || + NumTraits::IsComplex, + typename NumTraits::Real>::type + abs(const T& x) { + EIGEN_USING_STD_MATH(abs); + return abs(x); +} + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE + typename internal::enable_if::IsSigned || + NumTraits::IsComplex), + typename NumTraits::Real>::type + abs(const T& x) { + return x; +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_INTEGER_TYPES_UNARY(abs, abs) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(abs, fabs) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float abs(const float& x) { + return ::fabsf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double abs(const double& x) { + return ::fabs(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float abs(const std::complex& x) { + return ::hypotf(x.real(), x.imag()); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double abs( + const std::complex& x) { + return ::hypot(x.real(), x.imag()); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T exp(const T& x) { + EIGEN_USING_STD_MATH(exp); + return exp(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(exp, exp) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float exp(const float& x) { + return ::expf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double exp(const double& x) { + return ::exp(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex exp( + const std::complex& x) { + float com = ::expf(x.real()); + float res_real = com * ::cosf(x.imag()); + float res_imag = com * ::sinf(x.imag()); + return std::complex(res_real, res_imag); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex exp( + const std::complex& x) { + double com = ::exp(x.real()); + double res_real = com * ::cos(x.imag()); + double res_imag = com * ::sin(x.imag()); + return std::complex(res_real, res_imag); +} +#endif + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(expm1, Scalar) + expm1(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(expm1, Scalar)::run(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(expm1, expm1) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float expm1(const float& x) { + return ::expm1f(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double expm1(const double& x) { + return ::expm1(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T cos(const T& x) { + EIGEN_USING_STD_MATH(cos); + return cos(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(cos, cos) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float cos(const float& x) { + return ::cosf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double cos(const double& x) { + return ::cos(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T sin(const T& x) { + EIGEN_USING_STD_MATH(sin); + return sin(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(sin, sin) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float sin(const float& x) { + return ::sinf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double sin(const double& x) { + return ::sin(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T tan(const T& x) { + EIGEN_USING_STD_MATH(tan); + return tan(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(tan, tan) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float tan(const float& x) { + return ::tanf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double tan(const double& x) { + return ::tan(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T acos(const T& x) { + EIGEN_USING_STD_MATH(acos); + return acos(x); +} + +#if EIGEN_HAS_CXX11_MATH +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T acosh(const T& x) { + EIGEN_USING_STD_MATH(acosh); + return acosh(x); +} +#endif + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(acos, acos) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(acosh, acosh) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float acos(const float& x) { + return ::acosf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double acos(const double& x) { + return ::acos(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T asin(const T& x) { + EIGEN_USING_STD_MATH(asin); + return asin(x); +} + +#if EIGEN_HAS_CXX11_MATH +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T asinh(const T& x) { + EIGEN_USING_STD_MATH(asinh); + return asinh(x); +} +#endif + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(asin, asin) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(asinh, asinh) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float asin(const float& x) { + return ::asinf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double asin(const double& x) { + return ::asin(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T atan(const T& x) { + EIGEN_USING_STD_MATH(atan); + return atan(x); +} + +#if EIGEN_HAS_CXX11_MATH +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T atanh(const T& x) { + EIGEN_USING_STD_MATH(atanh); + return atanh(x); +} +#endif + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(atan, atan) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(atanh, atanh) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float atan(const float& x) { + return ::atanf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double atan(const double& x) { + return ::atan(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T cosh(const T& x) { + EIGEN_USING_STD_MATH(cosh); + return cosh(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(cosh, cosh) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float cosh(const float& x) { + return ::coshf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double cosh(const double& x) { + return ::cosh(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T sinh(const T& x) { + EIGEN_USING_STD_MATH(sinh); + return sinh(x); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(sinh, sinh) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float sinh(const float& x) { + return ::sinhf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double sinh(const double& x) { + return ::sinh(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T tanh(const T& x) { + EIGEN_USING_STD_MATH(tanh); + return tanh(x); +} + +#if (!defined(EIGEN_GPUCC)) && EIGEN_FAST_MATH && !defined(SYCL_DEVICE_ONLY) +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float tanh(float x) { + return internal::generic_fast_tanh_float(x); +} +#endif + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(tanh, tanh) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float tanh(const float& x) { + return ::tanhf(x); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double tanh(const double& x) { + return ::tanh(x); +} +#endif + +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T fmod(const T& a, const T& b) { + EIGEN_USING_STD_MATH(fmod); + return fmod(a, b); +} + +#if defined(SYCL_DEVICE_ONLY) +SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(fmod, fmod) +#endif + +#if defined(EIGEN_GPUCC) +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float fmod(const float& a, + const float& b) { + return ::fmodf(a, b); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double fmod(const double& a, + const double& b) { + return ::fmod(a, b); +} +#endif + +#if defined(SYCL_DEVICE_ONLY) +#undef SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY +#undef SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY +#undef SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_BINARY +#undef SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY +#undef SYCL_SPECIALIZE_INTEGER_TYPES_BINARY +#undef SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY +#undef SYCL_SPECIALIZE_FLOATING_TYPES_BINARY +#undef SYCL_SPECIALIZE_FLOATING_TYPES_UNARY +#undef SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE +#undef SYCL_SPECIALIZE_GEN_UNARY_FUNC +#undef SYCL_SPECIALIZE_UNARY_FUNC +#undef SYCL_SPECIALIZE_GEN1_BINARY_FUNC +#undef SYCL_SPECIALIZE_GEN2_BINARY_FUNC +#undef SYCL_SPECIALIZE_BINARY_FUNC +#endif + +} // end namespace numext + +namespace internal { + +template +EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex& x) { + return (numext::isfinite)(numext::real(x)) && + (numext::isfinite)(numext::imag(x)); +} + +template +EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex& x) { + return (numext::isnan)(numext::real(x)) || (numext::isnan)(numext::imag(x)); +} + +template +EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex& x) { + return ((numext::isinf)(numext::real(x)) || + (numext::isinf)(numext::imag(x))) && + (!(numext::isnan)(x)); +} + +/**************************************************************************** +* Implementation of fuzzy comparisons * +****************************************************************************/ + +template +struct scalar_fuzzy_default_impl {}; + +template +struct scalar_fuzzy_default_impl { + typedef typename NumTraits::Real RealScalar; + template + EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan( + const Scalar& x, const OtherScalar& y, const RealScalar& prec) { + return numext::abs(x) <= numext::abs(y) * prec; + } + EIGEN_DEVICE_FUNC + static inline bool isApprox(const Scalar& x, + const Scalar& y, + const RealScalar& prec) { + return numext::abs(x - y) <= + numext::mini(numext::abs(x), numext::abs(y)) * prec; + } + EIGEN_DEVICE_FUNC + static inline bool isApproxOrLessThan(const Scalar& x, + const Scalar& y, + const RealScalar& prec) { + return x <= y || isApprox(x, y, prec); + } +}; + +template +struct scalar_fuzzy_default_impl { + typedef typename NumTraits::Real RealScalar; + template + EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const Scalar& x, + const Scalar&, + const RealScalar&) { + return x == Scalar(0); + } + EIGEN_DEVICE_FUNC + static inline bool isApprox(const Scalar& x, + const Scalar& y, + const RealScalar&) { + return x == y; + } + EIGEN_DEVICE_FUNC + static inline bool isApproxOrLessThan(const Scalar& x, + const Scalar& y, + const RealScalar&) { + return x <= y; + } +}; + +template +struct scalar_fuzzy_default_impl { + typedef typename NumTraits::Real RealScalar; + template + EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan( + const Scalar& x, const OtherScalar& y, const RealScalar& prec) { + return numext::abs2(x) <= numext::abs2(y) * prec * prec; + } + EIGEN_DEVICE_FUNC + static inline bool isApprox(const Scalar& x, + const Scalar& y, + const RealScalar& prec) { + return numext::abs2(x - y) <= + numext::mini(numext::abs2(x), numext::abs2(y)) * prec * prec; + } +}; + +template +struct scalar_fuzzy_impl + : scalar_fuzzy_default_impl::IsComplex, + NumTraits::IsInteger> {}; + +template +EIGEN_DEVICE_FUNC inline bool isMuchSmallerThan( + const Scalar& x, + const OtherScalar& y, + const typename NumTraits::Real& precision = + NumTraits::dummy_precision()) { + return scalar_fuzzy_impl::template isMuchSmallerThan( + x, y, precision); +} + +template +EIGEN_DEVICE_FUNC inline bool isApprox( + const Scalar& x, + const Scalar& y, + const typename NumTraits::Real& precision = + NumTraits::dummy_precision()) { + return scalar_fuzzy_impl::isApprox(x, y, precision); +} + +template +EIGEN_DEVICE_FUNC inline bool isApproxOrLessThan( + const Scalar& x, + const Scalar& y, + const typename NumTraits::Real& precision = + NumTraits::dummy_precision()) { + return scalar_fuzzy_impl::isApproxOrLessThan(x, y, precision); +} + +/****************************************** +*** The special case of the bool type *** +******************************************/ + +template <> +struct random_impl { + static inline bool run() { return random(0, 1) == 0 ? false : true; } +}; + +template <> +struct scalar_fuzzy_impl { + typedef bool RealScalar; + + template + EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const bool& x, + const bool&, + const bool&) { + return !x; + } + + EIGEN_DEVICE_FUNC + static inline bool isApprox(bool x, bool y, bool) { return x == y; } + + EIGEN_DEVICE_FUNC + static inline bool isApproxOrLessThan(const bool& x, + const bool& y, + const bool&) { + return (!x) || y; + } +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_MATHFUNCTIONS_H