From e2f21465fea76a80966f12a20d0be36597f19b44 Mon Sep 17 00:00:00 2001 From: Antonio Sanchez Date: Wed, 2 Dec 2020 14:00:57 -0800 Subject: Special function implementations for half/bfloat16 packets. Current implementations fail to consider half-float packets, only half-float scalars. Added specializations for packets on AVX, AVX512 and NEON. Added tests to `special_packetmath`. The current `special_functions` tests would fail for half and bfloat16 due to lack of precision. The NEON tests also fail with precision issues and due to different handling of `sqrt(inf)`, so special functions bessel, ndtri have been disabled. Tested with AVX, AVX512. --- unsupported/Eigen/SpecialFunctions | 19 +- .../src/SpecialFunctions/BesselFunctionsImpl.h | 132 ++++---- .../SpecialFunctions/BesselFunctionsPacketMath.h | 36 +- .../src/SpecialFunctions/SpecialFunctionsImpl.h | 5 +- .../SpecialFunctions/arch/AVX/BesselFunctions.h | 46 +++ .../SpecialFunctions/arch/AVX/SpecialFunctions.h | 16 + .../SpecialFunctions/arch/AVX512/BesselFunctions.h | 46 +++ .../arch/AVX512/SpecialFunctions.h | 16 + .../arch/GPU/GpuSpecialFunctions.h | 369 --------------------- .../SpecialFunctions/arch/GPU/SpecialFunctions.h | 369 +++++++++++++++++++++ .../SpecialFunctions/arch/NEON/BesselFunctions.h | 54 +++ .../SpecialFunctions/arch/NEON/SpecialFunctions.h | 34 ++ unsupported/test/special_functions.cpp | 89 +++-- unsupported/test/special_packetmath.cpp | 57 ++-- 14 files changed, 767 insertions(+), 521 deletions(-) create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/AVX/BesselFunctions.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/AVX/SpecialFunctions.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/AVX512/BesselFunctions.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/AVX512/SpecialFunctions.h delete mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/GPU/SpecialFunctions.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/NEON/BesselFunctions.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/arch/NEON/SpecialFunctions.h (limited to 'unsupported') diff --git a/unsupported/Eigen/SpecialFunctions b/unsupported/Eigen/SpecialFunctions index dda6618de..f6a2460e6 100644 --- a/unsupported/Eigen/SpecialFunctions +++ b/unsupported/Eigen/SpecialFunctions @@ -61,23 +61,36 @@ namespace Eigen { } #include "src/SpecialFunctions/BesselFunctionsImpl.h" -#include "src/SpecialFunctions/BesselFunctionsPacketMath.h" #include "src/SpecialFunctions/BesselFunctionsBFloat16.h" #include "src/SpecialFunctions/BesselFunctionsHalf.h" +#include "src/SpecialFunctions/BesselFunctionsPacketMath.h" #include "src/SpecialFunctions/BesselFunctionsFunctors.h" #include "src/SpecialFunctions/BesselFunctionsArrayAPI.h" #include "src/SpecialFunctions/SpecialFunctionsImpl.h" #if defined(EIGEN_HIPCC) #include "src/SpecialFunctions/HipVectorCompatibility.h" #endif -#include "src/SpecialFunctions/SpecialFunctionsPacketMath.h" #include "src/SpecialFunctions/SpecialFunctionsBFloat16.h" #include "src/SpecialFunctions/SpecialFunctionsHalf.h" +#include "src/SpecialFunctions/SpecialFunctionsPacketMath.h" #include "src/SpecialFunctions/SpecialFunctionsFunctors.h" #include "src/SpecialFunctions/SpecialFunctionsArrayAPI.h" +#if defined EIGEN_VECTORIZE_AVX512 + #include "src/SpecialFunctions/arch/AVX/BesselFunctions.h" + #include "src/SpecialFunctions/arch/AVX/SpecialFunctions.h" + #include "src/SpecialFunctions/arch/AVX512/BesselFunctions.h" + #include "src/SpecialFunctions/arch/AVX512/SpecialFunctions.h" +#elif defined EIGEN_VECTORIZE_AVX + #include "src/SpecialFunctions/arch/AVX/BesselFunctions.h" + #include "src/SpecialFunctions/arch/AVX/SpecialFunctions.h" +#elif defined EIGEN_VECTORIZE_NEON + #include "src/SpecialFunctions/arch/NEON/BesselFunctions.h" + #include "src/SpecialFunctions/arch/NEON/SpecialFunctions.h" +#endif + #if defined EIGEN_VECTORIZE_GPU - #include "src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h" + #include "src/SpecialFunctions/arch/GPU/SpecialFunctions.h" #endif namespace Eigen { diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h index a9b6ad940..24812be1b 100644 --- a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h @@ -46,7 +46,7 @@ struct bessel_i0e_retval { typedef Scalar type; }; -template +template ::type> struct generic_i0e { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -201,11 +201,11 @@ struct generic_i0e { } }; -template +template struct bessel_i0e_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_i0e::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_i0e::run(x); } }; @@ -214,7 +214,7 @@ struct bessel_i0_retval { typedef Scalar type; }; -template +template ::type> struct generic_i0 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { @@ -224,11 +224,11 @@ struct generic_i0 { } }; -template +template struct bessel_i0_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_i0::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_i0::run(x); } }; @@ -237,7 +237,7 @@ struct bessel_i1e_retval { typedef Scalar type; }; -template +template ::type > struct generic_i1e { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -396,20 +396,20 @@ struct generic_i1e { } }; -template +template struct bessel_i1e_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_i1e::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_i1e::run(x); } }; -template +template struct bessel_i1_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_i1 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T& x) { @@ -419,20 +419,20 @@ struct generic_i1 { } }; -template +template struct bessel_i1_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_i1::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_i1::run(x); } }; -template +template struct bessel_k0e_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_k0e { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -582,20 +582,20 @@ struct generic_k0e { } }; -template +template struct bessel_k0e_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_k0e::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_k0e::run(x); } }; -template +template struct bessel_k0_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_k0 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -754,20 +754,20 @@ struct generic_k0 { } }; -template +template struct bessel_k0_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_k0::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_k0::run(x); } }; -template +template struct bessel_k1e_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_k1e { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -910,20 +910,20 @@ struct generic_k1e { } }; -template +template struct bessel_k1e_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_k1e::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_k1e::run(x); } }; -template +template struct bessel_k1_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_k1 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -1076,20 +1076,20 @@ struct generic_k1 { } }; -template +template struct bessel_k1_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_k1::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_k1::run(x); } }; -template +template struct bessel_j0_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_j0 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -1276,20 +1276,20 @@ struct generic_j0 { } }; -template +template struct bessel_j0_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_j0::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_j0::run(x); } }; -template +template struct bessel_y0_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_y0 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -1474,20 +1474,20 @@ struct generic_y0 { } }; -template +template struct bessel_y0_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_y0::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_y0::run(x); } }; -template +template struct bessel_j1_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_j1 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -1665,20 +1665,20 @@ struct generic_j1 { } }; -template +template struct bessel_j1_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_j1::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_j1::run(x); } }; -template +template struct bessel_y1_retval { - typedef Scalar type; + typedef T type; }; -template +template ::type> struct generic_y1 { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T run(const T&) { @@ -1868,11 +1868,11 @@ struct generic_y1 { } }; -template +template struct bessel_y1_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { - return generic_y1::run(x); + static EIGEN_STRONG_INLINE T run(const T x) { + return generic_y1::run(x); } }; diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h index efc6d9c8f..943d10f6a 100644 --- a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h +++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h @@ -19,8 +19,7 @@ namespace internal { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i0(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_i0; return generic_i0::run(x); + return numext::bessel_i0(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -28,8 +27,7 @@ Packet pbessel_i0(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i0e(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_i0e; return generic_i0e::run(x); + return numext::bessel_i0e(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -37,8 +35,7 @@ Packet pbessel_i0e(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i1(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_i1; return generic_i1::run(x); + return numext::bessel_i1(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -46,8 +43,7 @@ Packet pbessel_i1(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i1e(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_i1e; return generic_i1e::run(x); + return numext::bessel_i1e(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -55,8 +51,7 @@ Packet pbessel_i1e(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_j0(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_j0; return generic_j0::run(x); + return numext::bessel_j0(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -64,8 +59,7 @@ Packet pbessel_j0(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_j1(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_j1; return generic_j1::run(x); + return numext::bessel_j1(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -73,8 +67,7 @@ Packet pbessel_j1(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_y0(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_y0; return generic_y0::run(x); + return numext::bessel_y0(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -82,8 +75,7 @@ Packet pbessel_y0(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_y1(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_y1; return generic_y1::run(x); + return numext::bessel_y1(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -91,8 +83,7 @@ Packet pbessel_y1(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k0(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_k0; return generic_k0::run(x); + return numext::bessel_k0(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -100,8 +91,7 @@ Packet pbessel_k0(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k0e(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_k0e; return generic_k0e::run(x); + return numext::bessel_k0e(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -109,8 +99,7 @@ Packet pbessel_k0e(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k1(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_k1; return generic_k1::run(x); + return numext::bessel_k1(x); } /** \internal \returns the exponentially scaled modified Bessel function of @@ -118,8 +107,7 @@ Packet pbessel_k1(const Packet& x) { template EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k1e(const Packet& x) { - typedef typename unpacket_traits::type ScalarType; - using internal::generic_k1e; return generic_k1e::run(x); + return numext::bessel_k1e(x); } } // end namespace internal diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index 648eb053e..cfc13aff7 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -348,7 +348,7 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erf_float(const T& a_x) { template struct erf_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE T run(const T x) { + static EIGEN_STRONG_INLINE T run(const T& x) { return generic_fast_erf_float(x); } }; @@ -490,7 +490,8 @@ struct erfc_impl { template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T flipsign( const T& should_flipsign, const T& x) { - const T sign_mask = pset1(-0.0); + typedef typename unpacket_traits::type Scalar; + const T sign_mask = pset1(Scalar(-0.0)); T sign_bit = pand(should_flipsign, sign_mask); return pxor(sign_bit, x); } diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/AVX/BesselFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/AVX/BesselFunctions.h new file mode 100644 index 000000000..2d7669209 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/AVX/BesselFunctions.h @@ -0,0 +1,46 @@ +#ifndef EIGEN_AVX_BESSELFUNCTIONS_H +#define EIGEN_AVX_BESSELFUNCTIONS_H + +namespace Eigen { +namespace internal { + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_i0) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_i0) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_i0e) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_i0e) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_i1) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_i1) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_i1e) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_i1e) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_j0) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_j0) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_j1) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_j1) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_k0) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_k0) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_k0e) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_k0e) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_k1) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_k1) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_k1e) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_k1e) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_y0) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_y0) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pbessel_y1) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pbessel_y1) + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_AVX_BESSELFUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/AVX/SpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/AVX/SpecialFunctions.h new file mode 100644 index 000000000..35e62a8ac --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/AVX/SpecialFunctions.h @@ -0,0 +1,16 @@ +#ifndef EIGEN_AVX_SPECIALFUNCTIONS_H +#define EIGEN_AVX_SPECIALFUNCTIONS_H + +namespace Eigen { +namespace internal { + +F16_PACKET_FUNCTION(Packet8f, Packet8h, perf) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, perf) + +F16_PACKET_FUNCTION(Packet8f, Packet8h, pndtri) +BF16_PACKET_FUNCTION(Packet8f, Packet8bf, pndtri) + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_AVX_SPECIAL_FUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/AVX512/BesselFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/AVX512/BesselFunctions.h new file mode 100644 index 000000000..7dd3c3e5b --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/AVX512/BesselFunctions.h @@ -0,0 +1,46 @@ +#ifndef EIGEN_AVX512_BESSELFUNCTIONS_H +#define EIGEN_AVX512_BESSELFUNCTIONS_H + +namespace Eigen { +namespace internal { + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_i0) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_i0) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_i0e) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_i0e) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_i1) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_i1) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_i1e) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_i1e) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_j0) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_j0) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_j1) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_j1) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_k0) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_k0) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_k0e) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_k0e) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_k1) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_k1) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_k1e) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_k1e) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_y0) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_y0) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pbessel_y1) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pbessel_y1) + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_AVX512_BESSELFUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/AVX512/SpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/AVX512/SpecialFunctions.h new file mode 100644 index 000000000..79878f2b6 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/AVX512/SpecialFunctions.h @@ -0,0 +1,16 @@ +#ifndef EIGEN_AVX512_SPECIALFUNCTIONS_H +#define EIGEN_AVX512_SPECIALFUNCTIONS_H + +namespace Eigen { +namespace internal { + +F16_PACKET_FUNCTION(Packet16f, Packet16h, perf) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, perf) + +F16_PACKET_FUNCTION(Packet16f, Packet16h, pndtri) +BF16_PACKET_FUNCTION(Packet16f, Packet16bf, pndtri) + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_AVX512_SPECIAL_FUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h deleted file mode 100644 index dd3bf4dd1..000000000 --- a/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h +++ /dev/null @@ -1,369 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2014 Benoit Steiner -// -// 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_GPU_SPECIALFUNCTIONS_H -#define EIGEN_GPU_SPECIALFUNCTIONS_H - -namespace Eigen { - -namespace internal { - -// Make sure this is only available when targeting a GPU: we don't want to -// introduce conflicts between these packet_traits definitions and the ones -// we'll use on the host side (SSE, AVX, ...) -#if defined(EIGEN_GPUCC) && defined(EIGEN_USE_GPU) - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 plgamma(const float4& a) -{ - return make_float4(lgammaf(a.x), lgammaf(a.y), lgammaf(a.z), lgammaf(a.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 plgamma(const double2& a) -{ - using numext::lgamma; - return make_double2(lgamma(a.x), lgamma(a.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 pdigamma(const float4& a) -{ - using numext::digamma; - return make_float4(digamma(a.x), digamma(a.y), digamma(a.z), digamma(a.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 pdigamma(const double2& a) -{ - using numext::digamma; - return make_double2(digamma(a.x), digamma(a.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 pzeta(const float4& x, const float4& q) -{ - using numext::zeta; - return make_float4(zeta(x.x, q.x), zeta(x.y, q.y), zeta(x.z, q.z), zeta(x.w, q.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 pzeta(const double2& x, const double2& q) -{ - using numext::zeta; - return make_double2(zeta(x.x, q.x), zeta(x.y, q.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 ppolygamma(const float4& n, const float4& x) -{ - using numext::polygamma; - return make_float4(polygamma(n.x, x.x), polygamma(n.y, x.y), polygamma(n.z, x.z), polygamma(n.w, x.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 ppolygamma(const double2& n, const double2& x) -{ - using numext::polygamma; - return make_double2(polygamma(n.x, x.x), polygamma(n.y, x.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 perf(const float4& a) -{ - return make_float4(erff(a.x), erff(a.y), erff(a.z), erff(a.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 perf(const double2& a) -{ - using numext::erf; - return make_double2(erf(a.x), erf(a.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 perfc(const float4& a) -{ - using numext::erfc; - return make_float4(erfc(a.x), erfc(a.y), erfc(a.z), erfc(a.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 perfc(const double2& a) -{ - using numext::erfc; - return make_double2(erfc(a.x), erfc(a.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 pndtri(const float4& a) -{ - using numext::ndtri; - return make_float4(ndtri(a.x), ndtri(a.y), ndtri(a.z), ndtri(a.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 pndtri(const double2& a) -{ - using numext::ndtri; - return make_double2(ndtri(a.x), ndtri(a.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 pigamma(const float4& a, const float4& x) -{ - using numext::igamma; - return make_float4( - igamma(a.x, x.x), - igamma(a.y, x.y), - igamma(a.z, x.z), - igamma(a.w, x.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 pigamma(const double2& a, const double2& x) -{ - using numext::igamma; - return make_double2(igamma(a.x, x.x), igamma(a.y, x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pigamma_der_a( - const float4& a, const float4& x) { - using numext::igamma_der_a; - return make_float4(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y), - igamma_der_a(a.z, x.z), igamma_der_a(a.w, x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pigamma_der_a(const double2& a, const double2& x) { - using numext::igamma_der_a; - return make_double2(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pgamma_sample_der_alpha( - const float4& alpha, const float4& sample) { - using numext::gamma_sample_der_alpha; - return make_float4( - gamma_sample_der_alpha(alpha.x, sample.x), - gamma_sample_der_alpha(alpha.y, sample.y), - gamma_sample_der_alpha(alpha.z, sample.z), - gamma_sample_der_alpha(alpha.w, sample.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pgamma_sample_der_alpha(const double2& alpha, const double2& sample) { - using numext::gamma_sample_der_alpha; - return make_double2( - gamma_sample_der_alpha(alpha.x, sample.x), - gamma_sample_der_alpha(alpha.y, sample.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 pigammac(const float4& a, const float4& x) -{ - using numext::igammac; - return make_float4( - igammac(a.x, x.x), - igammac(a.y, x.y), - igammac(a.z, x.z), - igammac(a.w, x.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 pigammac(const double2& a, const double2& x) -{ - using numext::igammac; - return make_double2(igammac(a.x, x.x), igammac(a.y, x.y)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -float4 pbetainc(const float4& a, const float4& b, const float4& x) -{ - using numext::betainc; - return make_float4( - betainc(a.x, b.x, x.x), - betainc(a.y, b.y, x.y), - betainc(a.z, b.z, x.z), - betainc(a.w, b.w, x.w)); -} - -template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -double2 pbetainc(const double2& a, const double2& b, const double2& x) -{ - using numext::betainc; - return make_double2(betainc(a.x, b.x, x.x), betainc(a.y, b.y, x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i0e(const float4& x) { - using numext::bessel_i0e; - return make_float4(bessel_i0e(x.x), bessel_i0e(x.y), bessel_i0e(x.z), bessel_i0e(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_i0e(const double2& x) { - using numext::bessel_i0e; - return make_double2(bessel_i0e(x.x), bessel_i0e(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i0(const float4& x) { - using numext::bessel_i0; - return make_float4(bessel_i0(x.x), bessel_i0(x.y), bessel_i0(x.z), bessel_i0(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_i0(const double2& x) { - using numext::bessel_i0; - return make_double2(bessel_i0(x.x), bessel_i0(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i1e(const float4& x) { - using numext::bessel_i1e; - return make_float4(bessel_i1e(x.x), bessel_i1e(x.y), bessel_i1e(x.z), bessel_i1e(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_i1e(const double2& x) { - using numext::bessel_i1e; - return make_double2(bessel_i1e(x.x), bessel_i1e(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i1(const float4& x) { - using numext::bessel_i1; - return make_float4(bessel_i1(x.x), bessel_i1(x.y), bessel_i1(x.z), bessel_i1(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_i1(const double2& x) { - using numext::bessel_i1; - return make_double2(bessel_i1(x.x), bessel_i1(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k0e(const float4& x) { - using numext::bessel_k0e; - return make_float4(bessel_k0e(x.x), bessel_k0e(x.y), bessel_k0e(x.z), bessel_k0e(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_k0e(const double2& x) { - using numext::bessel_k0e; - return make_double2(bessel_k0e(x.x), bessel_k0e(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k0(const float4& x) { - using numext::bessel_k0; - return make_float4(bessel_k0(x.x), bessel_k0(x.y), bessel_k0(x.z), bessel_k0(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_k0(const double2& x) { - using numext::bessel_k0; - return make_double2(bessel_k0(x.x), bessel_k0(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k1e(const float4& x) { - using numext::bessel_k1e; - return make_float4(bessel_k1e(x.x), bessel_k1e(x.y), bessel_k1e(x.z), bessel_k1e(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_k1e(const double2& x) { - using numext::bessel_k1e; - return make_double2(bessel_k1e(x.x), bessel_k1e(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k1(const float4& x) { - using numext::bessel_k1; - return make_float4(bessel_k1(x.x), bessel_k1(x.y), bessel_k1(x.z), bessel_k1(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_k1(const double2& x) { - using numext::bessel_k1; - return make_double2(bessel_k1(x.x), bessel_k1(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_j0(const float4& x) { - using numext::bessel_j0; - return make_float4(bessel_j0(x.x), bessel_j0(x.y), bessel_j0(x.z), bessel_j0(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_j0(const double2& x) { - using numext::bessel_j0; - return make_double2(bessel_j0(x.x), bessel_j0(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_j1(const float4& x) { - using numext::bessel_j1; - return make_float4(bessel_j1(x.x), bessel_j1(x.y), bessel_j1(x.z), bessel_j1(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_j1(const double2& x) { - using numext::bessel_j1; - return make_double2(bessel_j1(x.x), bessel_j1(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_y0(const float4& x) { - using numext::bessel_y0; - return make_float4(bessel_y0(x.x), bessel_y0(x.y), bessel_y0(x.z), bessel_y0(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_y0(const double2& x) { - using numext::bessel_y0; - return make_double2(bessel_y0(x.x), bessel_y0(x.y)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_y1(const float4& x) { - using numext::bessel_y1; - return make_float4(bessel_y1(x.x), bessel_y1(x.y), bessel_y1(x.z), bessel_y1(x.w)); -} - -template <> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 -pbessel_y1(const double2& x) { - using numext::bessel_y1; - return make_double2(bessel_y1(x.x), bessel_y1(x.y)); -} - -#endif - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_GPU_SPECIALFUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/GPU/SpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/GPU/SpecialFunctions.h new file mode 100644 index 000000000..dd3bf4dd1 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/GPU/SpecialFunctions.h @@ -0,0 +1,369 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2014 Benoit Steiner +// +// 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_GPU_SPECIALFUNCTIONS_H +#define EIGEN_GPU_SPECIALFUNCTIONS_H + +namespace Eigen { + +namespace internal { + +// Make sure this is only available when targeting a GPU: we don't want to +// introduce conflicts between these packet_traits definitions and the ones +// we'll use on the host side (SSE, AVX, ...) +#if defined(EIGEN_GPUCC) && defined(EIGEN_USE_GPU) + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 plgamma(const float4& a) +{ + return make_float4(lgammaf(a.x), lgammaf(a.y), lgammaf(a.z), lgammaf(a.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 plgamma(const double2& a) +{ + using numext::lgamma; + return make_double2(lgamma(a.x), lgamma(a.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pdigamma(const float4& a) +{ + using numext::digamma; + return make_float4(digamma(a.x), digamma(a.y), digamma(a.z), digamma(a.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pdigamma(const double2& a) +{ + using numext::digamma; + return make_double2(digamma(a.x), digamma(a.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pzeta(const float4& x, const float4& q) +{ + using numext::zeta; + return make_float4(zeta(x.x, q.x), zeta(x.y, q.y), zeta(x.z, q.z), zeta(x.w, q.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pzeta(const double2& x, const double2& q) +{ + using numext::zeta; + return make_double2(zeta(x.x, q.x), zeta(x.y, q.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 ppolygamma(const float4& n, const float4& x) +{ + using numext::polygamma; + return make_float4(polygamma(n.x, x.x), polygamma(n.y, x.y), polygamma(n.z, x.z), polygamma(n.w, x.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 ppolygamma(const double2& n, const double2& x) +{ + using numext::polygamma; + return make_double2(polygamma(n.x, x.x), polygamma(n.y, x.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 perf(const float4& a) +{ + return make_float4(erff(a.x), erff(a.y), erff(a.z), erff(a.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 perf(const double2& a) +{ + using numext::erf; + return make_double2(erf(a.x), erf(a.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 perfc(const float4& a) +{ + using numext::erfc; + return make_float4(erfc(a.x), erfc(a.y), erfc(a.z), erfc(a.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 perfc(const double2& a) +{ + using numext::erfc; + return make_double2(erfc(a.x), erfc(a.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pndtri(const float4& a) +{ + using numext::ndtri; + return make_float4(ndtri(a.x), ndtri(a.y), ndtri(a.z), ndtri(a.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pndtri(const double2& a) +{ + using numext::ndtri; + return make_double2(ndtri(a.x), ndtri(a.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pigamma(const float4& a, const float4& x) +{ + using numext::igamma; + return make_float4( + igamma(a.x, x.x), + igamma(a.y, x.y), + igamma(a.z, x.z), + igamma(a.w, x.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pigamma(const double2& a, const double2& x) +{ + using numext::igamma; + return make_double2(igamma(a.x, x.x), igamma(a.y, x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pigamma_der_a( + const float4& a, const float4& x) { + using numext::igamma_der_a; + return make_float4(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y), + igamma_der_a(a.z, x.z), igamma_der_a(a.w, x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pigamma_der_a(const double2& a, const double2& x) { + using numext::igamma_der_a; + return make_double2(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pgamma_sample_der_alpha( + const float4& alpha, const float4& sample) { + using numext::gamma_sample_der_alpha; + return make_float4( + gamma_sample_der_alpha(alpha.x, sample.x), + gamma_sample_der_alpha(alpha.y, sample.y), + gamma_sample_der_alpha(alpha.z, sample.z), + gamma_sample_der_alpha(alpha.w, sample.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pgamma_sample_der_alpha(const double2& alpha, const double2& sample) { + using numext::gamma_sample_der_alpha; + return make_double2( + gamma_sample_der_alpha(alpha.x, sample.x), + gamma_sample_der_alpha(alpha.y, sample.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pigammac(const float4& a, const float4& x) +{ + using numext::igammac; + return make_float4( + igammac(a.x, x.x), + igammac(a.y, x.y), + igammac(a.z, x.z), + igammac(a.w, x.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pigammac(const double2& a, const double2& x) +{ + using numext::igammac; + return make_double2(igammac(a.x, x.x), igammac(a.y, x.y)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +float4 pbetainc(const float4& a, const float4& b, const float4& x) +{ + using numext::betainc; + return make_float4( + betainc(a.x, b.x, x.x), + betainc(a.y, b.y, x.y), + betainc(a.z, b.z, x.z), + betainc(a.w, b.w, x.w)); +} + +template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE +double2 pbetainc(const double2& a, const double2& b, const double2& x) +{ + using numext::betainc; + return make_double2(betainc(a.x, b.x, x.x), betainc(a.y, b.y, x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i0e(const float4& x) { + using numext::bessel_i0e; + return make_float4(bessel_i0e(x.x), bessel_i0e(x.y), bessel_i0e(x.z), bessel_i0e(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_i0e(const double2& x) { + using numext::bessel_i0e; + return make_double2(bessel_i0e(x.x), bessel_i0e(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i0(const float4& x) { + using numext::bessel_i0; + return make_float4(bessel_i0(x.x), bessel_i0(x.y), bessel_i0(x.z), bessel_i0(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_i0(const double2& x) { + using numext::bessel_i0; + return make_double2(bessel_i0(x.x), bessel_i0(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i1e(const float4& x) { + using numext::bessel_i1e; + return make_float4(bessel_i1e(x.x), bessel_i1e(x.y), bessel_i1e(x.z), bessel_i1e(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_i1e(const double2& x) { + using numext::bessel_i1e; + return make_double2(bessel_i1e(x.x), bessel_i1e(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_i1(const float4& x) { + using numext::bessel_i1; + return make_float4(bessel_i1(x.x), bessel_i1(x.y), bessel_i1(x.z), bessel_i1(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_i1(const double2& x) { + using numext::bessel_i1; + return make_double2(bessel_i1(x.x), bessel_i1(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k0e(const float4& x) { + using numext::bessel_k0e; + return make_float4(bessel_k0e(x.x), bessel_k0e(x.y), bessel_k0e(x.z), bessel_k0e(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_k0e(const double2& x) { + using numext::bessel_k0e; + return make_double2(bessel_k0e(x.x), bessel_k0e(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k0(const float4& x) { + using numext::bessel_k0; + return make_float4(bessel_k0(x.x), bessel_k0(x.y), bessel_k0(x.z), bessel_k0(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_k0(const double2& x) { + using numext::bessel_k0; + return make_double2(bessel_k0(x.x), bessel_k0(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k1e(const float4& x) { + using numext::bessel_k1e; + return make_float4(bessel_k1e(x.x), bessel_k1e(x.y), bessel_k1e(x.z), bessel_k1e(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_k1e(const double2& x) { + using numext::bessel_k1e; + return make_double2(bessel_k1e(x.x), bessel_k1e(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_k1(const float4& x) { + using numext::bessel_k1; + return make_float4(bessel_k1(x.x), bessel_k1(x.y), bessel_k1(x.z), bessel_k1(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_k1(const double2& x) { + using numext::bessel_k1; + return make_double2(bessel_k1(x.x), bessel_k1(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_j0(const float4& x) { + using numext::bessel_j0; + return make_float4(bessel_j0(x.x), bessel_j0(x.y), bessel_j0(x.z), bessel_j0(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_j0(const double2& x) { + using numext::bessel_j0; + return make_double2(bessel_j0(x.x), bessel_j0(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_j1(const float4& x) { + using numext::bessel_j1; + return make_float4(bessel_j1(x.x), bessel_j1(x.y), bessel_j1(x.z), bessel_j1(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_j1(const double2& x) { + using numext::bessel_j1; + return make_double2(bessel_j1(x.x), bessel_j1(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_y0(const float4& x) { + using numext::bessel_y0; + return make_float4(bessel_y0(x.x), bessel_y0(x.y), bessel_y0(x.z), bessel_y0(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_y0(const double2& x) { + using numext::bessel_y0; + return make_double2(bessel_y0(x.x), bessel_y0(x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pbessel_y1(const float4& x) { + using numext::bessel_y1; + return make_float4(bessel_y1(x.x), bessel_y1(x.y), bessel_y1(x.z), bessel_y1(x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pbessel_y1(const double2& x) { + using numext::bessel_y1; + return make_double2(bessel_y1(x.x), bessel_y1(x.y)); +} + +#endif + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_GPU_SPECIALFUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/NEON/BesselFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/NEON/BesselFunctions.h new file mode 100644 index 000000000..67433b057 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/NEON/BesselFunctions.h @@ -0,0 +1,54 @@ +#ifndef EIGEN_NEON_BESSELFUNCTIONS_H +#define EIGEN_NEON_BESSELFUNCTIONS_H + +namespace Eigen { +namespace internal { + +#if EIGEN_HAS_ARM64_FP16_VECTOR_ARITHMETIC + +#define NEON_HALF_TO_FLOAT_FUNCTIONS(METHOD) \ +template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \ +Packet8hf METHOD(const Packet8hf& x) { \ + const Packet4f lo = METHOD(vcvt_f32_f16(vget_low_f16(x))); \ + const Packet4f hi = METHOD(vcvt_f32_f16(vget_high_f16(x))); \ + return vcombine_f16(vcvt_f16_f32(lo), vcvt_f16_f32(hi)); \ +} \ + \ +template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \ +Packet4hf METHOD(const Packet4hf& x) { \ + return vcvt_f16_f32(METHOD(vcvt_f32_f16(x))); \ +} + +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_i0) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_i0e) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_i1) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_i1e) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_j0) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_j1) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_k0) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_k0e) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_k1) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_k1e) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_y0) +NEON_HALF_TO_FLOAT_FUNCTIONS(pbessel_y1) + +#undef NEON_HALF_TO_FLOAT_FUNCTIONS +#endif + +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_i0) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_i0e) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_i1) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_i1e) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_j0) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_j1) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_k0) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_k0e) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_k1) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_k1e) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_y0) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pbessel_y1) + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_NEON_BESSELFUNCTIONS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/NEON/SpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/NEON/SpecialFunctions.h new file mode 100644 index 000000000..f8dda28fc --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/arch/NEON/SpecialFunctions.h @@ -0,0 +1,34 @@ +#ifndef EIGEN_NEON_SPECIALFUNCTIONS_H +#define EIGEN_NEON_SPECIALFUNCTIONS_H + +namespace Eigen { +namespace internal { + +#ifdef EIGEN_HAS_ARM64_FP16_VECTOR_ARITHMETIC + +#define NEON_HALF_TO_FLOAT_FUNCTIONS(METHOD) \ +template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \ +Packet8hf METHOD(const Packet8hf& x) { \ + const Packet4f lo = METHOD(vcvt_f32_f16(vget_low_f16(x))); \ + const Packet4f hi = METHOD(vcvt_f32_f16(vget_high_f16(x))); \ + return vcombine_f16(vcvt_f16_f32(lo), vcvt_f16_f32(hi)); \ +} \ + \ +template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \ +Packet4hf METHOD(const Packet4hf& x) { \ + return vcvt_f16_f32(METHOD(vcvt_f32_f16(x))); \ +} + +NEON_HALF_TO_FLOAT_FUNCTIONS(perf) +NEON_HALF_TO_FLOAT_FUNCTIONS(pndtri) + +#undef NEON_HALF_TO_FLOAT_FUNCTIONS +#endif + +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, perf) +BF16_PACKET_FUNCTION(Packet4f, Packet4bf, pndtri) + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_NEON_SPECIALFUNCTIONS_H diff --git a/unsupported/test/special_functions.cpp b/unsupported/test/special_functions.cpp index 1027272d5..56848daa5 100644 --- a/unsupported/test/special_functions.cpp +++ b/unsupported/test/special_functions.cpp @@ -11,6 +11,17 @@ #include "main.h" #include "../Eigen/SpecialFunctions" +// Hack to allow "implicit" conversions from double to Scalar via comma-initialization. +template +Eigen::CommaInitializer operator<<(Eigen::DenseBase& dense, double v) { + return (dense << static_cast(v)); +} + +template +Eigen::CommaInitializer& operator,(Eigen::CommaInitializer& ci, double v) { + return (ci, static_cast(v)); +} + template void verify_component_wise(const X& x, const Y& y) { @@ -65,8 +76,8 @@ template void array_special_functions() // igamma(a, x) = gamma(a, x) / Gamma(a) // where Gamma and gamma are considered the standard unnormalized // upper and lower incomplete gamma functions, respectively. - ArrayType a = m1.abs() + 2; - ArrayType x = m2.abs() + 2; + ArrayType a = m1.abs() + Scalar(2); + ArrayType x = m2.abs() + Scalar(2); ArrayType zero = ArrayType::Zero(rows, cols); ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0)); ArrayType a_m1 = a - one; @@ -83,18 +94,18 @@ template void array_special_functions() VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp()); // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) - VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).exp()); + VERIFY_IS_APPROX(Gamma_a_x, (a - Scalar(1)) * Gamma_a_m1_x + x.pow(a-Scalar(1)) * (-x).exp()); // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) - VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).exp()); + VERIFY_IS_APPROX(gamma_a_x, (a - Scalar(1)) * gamma_a_m1_x - x.pow(a-Scalar(1)) * (-x).exp()); } { // Verify for large a and x that values are between 0 and 1. ArrayType m1 = ArrayType::Random(rows,cols); ArrayType m2 = ArrayType::Random(rows,cols); - Scalar max_exponent = std::numeric_limits::max_exponent10; - ArrayType a = m1.abs() * pow(10., max_exponent - 1); - ArrayType x = m2.abs() * pow(10., max_exponent - 1); + int max_exponent = std::numeric_limits::max_exponent10; + ArrayType a = m1.abs() * Scalar(pow(10., max_exponent - 1)); + ArrayType x = m2.abs() * Scalar(pow(10., max_exponent - 1)); for (int i = 0; i < a.size(); ++i) { Scalar igam = numext::igamma(a(i), x(i)); VERIFY(0 <= igam); @@ -108,27 +119,37 @@ template void array_special_functions() Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; // location i*6+j corresponds to a_s[i], x_s[j]. - Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan}, - {0.0, 0.6321205588285578, 0.7768698398515702, - 0.9816843611112658, 9.999500016666262e-05, 1.0}, - {0.0, 0.4275932955291202, 0.608374823728911, - 0.9539882943107686, 7.522076445089201e-07, 1.0}, - {0.0, 0.01898815687615381, 0.06564245437845008, - 0.5665298796332909, 4.166333347221828e-18, 1.0}, - {0.0, 0.9999780593618628, 0.9999899967080838, - 0.9999996219837988, 0.9991370418689945, 1.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}}; - Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan}, - {1.0, 0.36787944117144233, 0.22313016014842982, - 0.018315638888734182, 0.9999000049998333, 0.0}, - {1.0, 0.5724067044708798, 0.3916251762710878, - 0.04601170568923136, 0.9999992477923555, 0.0}, - {1.0, 0.9810118431238462, 0.9343575456215499, - 0.4334701203667089, 1.0, 0.0}, - {1.0, 2.1940638138146658e-05, 1.0003291916285e-05, - 3.7801620118431334e-07, 0.0008629581310054535, - 0.0}, - {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}}; + Scalar igamma_s[][6] = { + {Scalar(0.0), nan, nan, nan, nan, nan}, + {Scalar(0.0), Scalar(0.6321205588285578), Scalar(0.7768698398515702), + Scalar(0.9816843611112658), Scalar(9.999500016666262e-05), + Scalar(1.0)}, + {Scalar(0.0), Scalar(0.4275932955291202), Scalar(0.608374823728911), + Scalar(0.9539882943107686), Scalar(7.522076445089201e-07), + Scalar(1.0)}, + {Scalar(0.0), Scalar(0.01898815687615381), + Scalar(0.06564245437845008), Scalar(0.5665298796332909), + Scalar(4.166333347221828e-18), Scalar(1.0)}, + {Scalar(0.0), Scalar(0.9999780593618628), Scalar(0.9999899967080838), + Scalar(0.9999996219837988), Scalar(0.9991370418689945), Scalar(1.0)}, + {Scalar(0.0), Scalar(0.0), Scalar(0.0), Scalar(0.0), Scalar(0.0), + Scalar(0.5042041932513908)}}; + Scalar igammac_s[][6] = { + {nan, nan, nan, nan, nan, nan}, + {Scalar(1.0), Scalar(0.36787944117144233), + Scalar(0.22313016014842982), Scalar(0.018315638888734182), + Scalar(0.9999000049998333), Scalar(0.0)}, + {Scalar(1.0), Scalar(0.5724067044708798), Scalar(0.3916251762710878), + Scalar(0.04601170568923136), Scalar(0.9999992477923555), + Scalar(0.0)}, + {Scalar(1.0), Scalar(0.9810118431238462), Scalar(0.9343575456215499), + Scalar(0.4334701203667089), Scalar(1.0), Scalar(0.0)}, + {Scalar(1.0), Scalar(2.1940638138146658e-05), + Scalar(1.0003291916285e-05), Scalar(3.7801620118431334e-07), + Scalar(0.0008629581310054535), Scalar(0.0)}, + {Scalar(1.0), Scalar(1.0), Scalar(1.0), Scalar(1.0), Scalar(1.0), + Scalar(0.49579580674813944)}}; + for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) { if ((std::isnan)(igamma_s[i][j])) { @@ -162,8 +183,8 @@ template void array_special_functions() ArrayType m1 = ArrayType::Random(32); using std::sqrt; - ArrayType cdf_val = (m1 / sqrt(2.)).erf(); - cdf_val = (cdf_val + 1.) / 2.; + ArrayType cdf_val = (m1 / Scalar(sqrt(2.))).erf(); + cdf_val = (cdf_val + Scalar(1)) / Scalar(2); verify_component_wise(cdf_val.ndtri(), m1);); } @@ -190,7 +211,6 @@ template void array_special_functions() CALL_SUBTEST( res = digamma(x); verify_component_wise(res, ref); ); } - #if EIGEN_HAS_C99_MATH { ArrayType n(11), x(11), res(11), ref(11); @@ -323,8 +343,8 @@ template void array_special_functions() ArrayType m3 = ArrayType::Random(32); ArrayType one = ArrayType::Constant(32, Scalar(1.0)); const Scalar eps = std::numeric_limits::epsilon(); - ArrayType a = (m1 * 4.0).exp(); - ArrayType b = (m2 * 4.0).exp(); + ArrayType a = (m1 * Scalar(4)).exp(); + ArrayType b = (m2 * Scalar(4)).exp(); ArrayType x = m3.abs(); // betainc(a, 1, x) == x**a @@ -471,4 +491,7 @@ EIGEN_DECLARE_TEST(special_functions) { CALL_SUBTEST_1(array_special_functions()); CALL_SUBTEST_2(array_special_functions()); + // TODO(cantonios): half/bfloat16 don't have enough precision to reproduce results above. + // CALL_SUBTEST_3(array_special_functions>()); + // CALL_SUBTEST_4(array_special_functions>()); } diff --git a/unsupported/test/special_packetmath.cpp b/unsupported/test/special_packetmath.cpp index 87b87351b..4f426eb8a 100644 --- a/unsupported/test/special_packetmath.cpp +++ b/unsupported/test/special_packetmath.cpp @@ -8,6 +8,7 @@ // 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/. +#include #include "packetmath_test_shared.h" #include "../Eigen/SpecialFunctions" @@ -43,42 +44,48 @@ template void packetmath_real() } { for (int i=0; i(0,1); + data1[i] = internal::random(Scalar(0),Scalar(1)); } CHECK_CWISE1_IF(internal::packet_traits::HasNdtri, numext::ndtri, internal::pndtri); } #endif // EIGEN_HAS_C99_MATH // For bessel_i*e and bessel_j*, the valid range is negative reals. - for (int i=0; i(-1,1) * std::pow(Scalar(10), internal::random(-6,6)); - data2[i] = internal::random(-1,1) * std::pow(Scalar(10), internal::random(-6,6)); - } + const int max_exponent = numext::mini(std::numeric_limits::max_exponent10-1, 6); + for (int i=0; i(Scalar(-1),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random(Scalar(-max_exponent),Scalar(max_exponent)))); + data2[i] = internal::random(Scalar(-1),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random(Scalar(-max_exponent),Scalar(max_exponent)))); + } - CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e); - CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e); - CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0); - CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1); + CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e); + CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e); + CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0); + CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1); + } // Use a smaller data range for the bessel_i* as these can become very large. // Following #1693, we also restrict this range further to avoid inf's due to // differences in pexp and exp. for (int i=0; i(0.01,1) * std::pow( - Scalar(9), internal::random(-1,2)); - data2[i] = internal::random(0.01,1) * std::pow( - Scalar(9), internal::random(-1,2)); + data1[i] = internal::random(Scalar(0.01),Scalar(1)) * + Scalar(std::pow(Scalar(9), internal::random(Scalar(-1),Scalar(2)))); + data2[i] = internal::random(Scalar(0.01),Scalar(1)) * + Scalar(std::pow(Scalar(9), internal::random(Scalar(-1),Scalar(2)))); } CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0); CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1); // y_i, and k_i are valid for x > 0. - for (int i=0; i(0.01,1) * std::pow(Scalar(10), internal::random(-2,5)); - data2[i] = internal::random(0.01,1) * std::pow(Scalar(10), internal::random(-2,5)); + const int max_exponent = numext::mini(std::numeric_limits::max_exponent10-1, 5); + for (int i=0; i(Scalar(0.01),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random(Scalar(-2),Scalar(max_exponent)))); + data2[i] = internal::random(Scalar(0.01),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random(Scalar(-2),Scalar(max_exponent)))); + } } // TODO(srvasude): Re-enable this test once properly investigated why the @@ -91,20 +98,20 @@ template void packetmath_real() // Following #1693, we restrict the range for exp to avoid zeroing out too // fast. for (int i=0; i(0.01,1) * std::pow( - Scalar(9), internal::random(-1,2)); - data2[i] = internal::random(0.01,1) * std::pow( - Scalar(9), internal::random(-1,2)); + data1[i] = internal::random(Scalar(0.01),Scalar(1)) * + Scalar(std::pow(Scalar(9), internal::random(Scalar(-1),Scalar(2)))); + data2[i] = internal::random(Scalar(0.01),Scalar(1)) * + Scalar(std::pow(Scalar(9), internal::random(Scalar(-1),Scalar(2)))); } CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0, internal::pbessel_k0); CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1, internal::pbessel_k1); for (int i=0; i(0.01,1) * std::pow( - Scalar(10), internal::random(-1,2)); - data2[i] = internal::random(0.01,1) * std::pow( - Scalar(10), internal::random(-1,2)); + data1[i] = internal::random(Scalar(0.01),Scalar(1)) * + Scalar(std::pow(Scalar(10), internal::random(Scalar(-1),Scalar(2)))); + data2[i] = internal::random(Scalar(0.01),Scalar(1)) * + Scalar(std::pow(Scalar(10), internal::random(Scalar(-1),Scalar(2)))); } #if EIGEN_HAS_C99_MATH && (__cplusplus > 199711L) @@ -135,6 +142,8 @@ EIGEN_DECLARE_TEST(special_packetmath) CALL_SUBTEST_1( test::runner::run() ); CALL_SUBTEST_2( test::runner::run() ); + CALL_SUBTEST_3( test::runner::run() ); + CALL_SUBTEST_4( test::runner::run() ); g_first_pass = false; } } -- cgit v1.2.3