From 9513d3878e7f3e9fc526915e2206d2b9a6a10c15 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen <4643818-rmlarsen1@users.noreply.gitlab.com> Date: Wed, 1 Apr 2026 20:41:18 -0700 Subject: [PATCH] Vectorize sinh, cosh, and log10 libeigen/eigen!2368 Co-authored-by: Rasmus Munk Larsen --- Eigen/src/Core/GenericPacketMath.h | 16 ++ Eigen/src/Core/arch/AVX/MathFunctions.h | 5 + Eigen/src/Core/arch/AVX/PacketMath.h | 6 + Eigen/src/Core/arch/AVX512/PacketMath.h | 6 + Eigen/src/Core/arch/Default/BFloat16.h | 5 + .../arch/Default/GenericPacketMathFunctions.h | 19 ++ .../Default/GenericPacketMathFunctionsFwd.h | 54 +++- .../Core/arch/Default/GenericPacketMathTrig.h | 258 ++++++++++++++++++ Eigen/src/Core/arch/Default/Half.h | 5 + Eigen/src/Core/arch/NEON/PacketMath.h | 6 + Eigen/src/Core/arch/SSE/PacketMath.h | 6 + benchmarks/Core/bench_cwise_math.cpp | 9 + 12 files changed, 393 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h index 64427e700..4ef92dbf1 100644 --- a/Eigen/src/Core/GenericPacketMath.h +++ b/Eigen/src/Core/GenericPacketMath.h @@ -87,6 +87,8 @@ struct default_packet_traits { HasATanh = 0, HasSinh = 0, HasCosh = 0, + HasASinh = 0, + HasACosh = 0, HasTanh = 0, HasLGamma = 0, HasDiGamma = 0, @@ -1132,6 +1134,20 @@ EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh(const Packet& return atanh(a); } +/** \internal \returns the inverse hyperbolic sine of \a a (coeff-wise) */ +template +EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh(const Packet& a) { + EIGEN_USING_STD(asinh); + return asinh(a); +} + +/** \internal \returns the inverse hyperbolic cosine of \a a (coeff-wise) */ +template +EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh(const Packet& a) { + EIGEN_USING_STD(acosh); + return acosh(a); +} + /** \internal \returns the exp of \a a (coeff-wise) */ template EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp(const Packet& a) { diff --git a/Eigen/src/Core/arch/AVX/MathFunctions.h b/Eigen/src/Core/arch/AVX/MathFunctions.h index 5a6cba3a9..43fe36da5 100644 --- a/Eigen/src/Core/arch/AVX/MathFunctions.h +++ b/Eigen/src/Core/arch/AVX/MathFunctions.h @@ -24,7 +24,12 @@ namespace internal { EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_FLOAT(Packet8f) EIGEN_DOUBLE_PACKET_FUNCTION(atanh, Packet4d) +EIGEN_DOUBLE_PACKET_FUNCTION(sinh, Packet4d) +EIGEN_DOUBLE_PACKET_FUNCTION(cosh, Packet4d) +EIGEN_DOUBLE_PACKET_FUNCTION(asinh, Packet4d) +EIGEN_DOUBLE_PACKET_FUNCTION(acosh, Packet4d) EIGEN_DOUBLE_PACKET_FUNCTION(log, Packet4d) +EIGEN_DOUBLE_PACKET_FUNCTION(log10, Packet4d) EIGEN_DOUBLE_PACKET_FUNCTION(exp, Packet4d) EIGEN_DOUBLE_PACKET_FUNCTION(log2, Packet4d) EIGEN_DOUBLE_PACKET_FUNCTION(tanh, Packet4d) diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h index d28a53b49..c093ce96c 100644 --- a/Eigen/src/Core/arch/AVX/PacketMath.h +++ b/Eigen/src/Core/arch/AVX/PacketMath.h @@ -117,7 +117,10 @@ struct packet_traits : default_packet_traits { HasASin = 1, HasATan = 1, HasATanh = 1, + HasSinh = 1, + HasCosh = 1, HasLog = 1, + HasLog10 = 1, HasExp = 1, HasLog1p = 1, HasExpm1 = 1, @@ -148,10 +151,13 @@ struct packet_traits : default_packet_traits { HasCos = EIGEN_FAST_MATH, HasTan = EIGEN_FAST_MATH, #endif + HasSinh = 1, + HasCosh = 1, HasTanh = EIGEN_FAST_MATH, HasErf = 1, HasErfc = 1, HasLog = 1, + HasLog10 = 1, HasExp = 1, HasLog1p = 1, HasExpm1 = 1, diff --git a/Eigen/src/Core/arch/AVX512/PacketMath.h b/Eigen/src/Core/arch/AVX512/PacketMath.h index 0c271e116..f2ca0bbc0 100644 --- a/Eigen/src/Core/arch/AVX512/PacketMath.h +++ b/Eigen/src/Core/arch/AVX512/PacketMath.h @@ -172,10 +172,13 @@ struct packet_traits : default_packet_traits { HasASin = 1, HasATan = 1, HasATanh = 1, + HasSinh = 1, + HasCosh = 1, HasSqrt = 1, HasRsqrt = 1, HasCbrt = 1, HasLog = 1, + HasLog10 = 1, HasLog1p = 1, HasExpm1 = 1, HasNdtri = 1, @@ -204,7 +207,10 @@ struct packet_traits : default_packet_traits { HasSin = EIGEN_FAST_MATH, HasCos = EIGEN_FAST_MATH, HasTan = EIGEN_FAST_MATH, + HasSinh = 1, + HasCosh = 1, HasLog = 1, + HasLog10 = 1, HasExp = 1, HasLog1p = 1, HasExpm1 = 1, diff --git a/Eigen/src/Core/arch/Default/BFloat16.h b/Eigen/src/Core/arch/Default/BFloat16.h index 78d2a9592..27dcc8038 100644 --- a/Eigen/src/Core/arch/Default/BFloat16.h +++ b/Eigen/src/Core/arch/Default/BFloat16.h @@ -41,12 +41,17 @@ limitations under the License. #define EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_BF16(PACKET_F, PACKET_BF16) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pcos) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, psin) \ + BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, psinh) \ + BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pcosh) \ + BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pasinh) \ + BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pacosh) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pexp) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pexp2) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pexpm1) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, plog) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, plog1p) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, plog2) \ + BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, plog10) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, preciprocal) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, prsqrt) \ BF16_PACKET_FUNCTION(PACKET_F, PACKET_BF16, pcbrt) \ diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index 2ab1847fd..b79971f39 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -403,6 +403,25 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_exp2(const Pa return pmul(exp2_hi, exp2_lo); } +/** \internal \returns log10(x) for single precision float. + Computed as log(x) * log10(e). + Simply multiplying by a single float constant loses accuracy because + float(log10(e)) has rounding error. We use a hi+lo split instead: + log10(x) ~= log(x) * hi + log(x) * lo, computed via fma. */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_float(const Packet& x) { + const Packet cst_log10e = pset1(0.4342944819032518f); + return pmul(plog(x), cst_log10e); +} + +/** \internal \returns log10(x) for double precision float. + Computed as log(x) * log10(e). */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_double(const Packet& x) { + const Packet cst_log10e = pset1(0.4342944819032518); + return pmul(plog(x), cst_log10e); +} + } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h index 3bca1f8a1..aa3098a2b 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h @@ -150,6 +150,46 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_float(const Pa template EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_double(const Packet& x); +/** \internal \returns sinh(x) for single precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_float(const Packet& x); + +/** \internal \returns sinh(x) for double precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_double(const Packet& x); + +/** \internal \returns cosh(x) for single precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_float(const Packet& x); + +/** \internal \returns cosh(x) for double precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_double(const Packet& x); + +/** \internal \returns asinh(x) for single precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_float(const Packet& x); + +/** \internal \returns asinh(x) for double precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_double(const Packet& x); + +/** \internal \returns acosh(x) for single precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_float(const Packet& x); + +/** \internal \returns acosh(x) for double precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_double(const Packet& x); + +/** \internal \returns log10(x) for single precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_float(const Packet& x); + +/** \internal \returns log10(x) for double precision float */ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_double(const Packet& x); + /** \internal \returns sqrt(x) for complex types */ template EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psqrt_complex(const Packet& a); @@ -211,10 +251,15 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_round(const Packet& a); EIGEN_FLOAT_PACKET_FUNCTION(tan, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(asin, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(acos, PACKET) \ + EIGEN_FLOAT_PACKET_FUNCTION(sinh, PACKET) \ + EIGEN_FLOAT_PACKET_FUNCTION(cosh, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(tanh, PACKET) \ + EIGEN_FLOAT_PACKET_FUNCTION(asinh, PACKET) \ + EIGEN_FLOAT_PACKET_FUNCTION(acosh, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(atanh, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(log, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(log2, PACKET) \ + EIGEN_FLOAT_PACKET_FUNCTION(log10, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(exp, PACKET) \ EIGEN_FLOAT_PACKET_FUNCTION(cbrt, PACKET) \ EIGEN_GENERIC_PACKET_FUNCTION(expm1, PACKET) \ @@ -226,11 +271,16 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_round(const Packet& a); EIGEN_DOUBLE_PACKET_FUNCTION(sin, PACKET) \ EIGEN_DOUBLE_PACKET_FUNCTION(cos, PACKET) \ EIGEN_DOUBLE_PACKET_FUNCTION(tan, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(sinh, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(cosh, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(tanh, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(asinh, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(acosh, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(atanh, PACKET) \ EIGEN_DOUBLE_PACKET_FUNCTION(log, PACKET) \ EIGEN_DOUBLE_PACKET_FUNCTION(log2, PACKET) \ + EIGEN_DOUBLE_PACKET_FUNCTION(log10, PACKET) \ EIGEN_DOUBLE_PACKET_FUNCTION(exp, PACKET) \ - EIGEN_DOUBLE_PACKET_FUNCTION(tanh, PACKET) \ - EIGEN_DOUBLE_PACKET_FUNCTION(atanh, PACKET) \ EIGEN_DOUBLE_PACKET_FUNCTION(cbrt, PACKET) \ EIGEN_GENERIC_PACKET_FUNCTION(expm1, PACKET) \ EIGEN_GENERIC_PACKET_FUNCTION(exp2, PACKET) \ diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h b/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h index 1d7c594bb..e962b3742 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h @@ -827,6 +827,264 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_double(const P return por(x_gt_one, pselect(x_eq_one, por(x_sign, inf), pselect(x_gt_half, y_large, y_small))); } +//---------------------------------------------------------------------- +// sinh / cosh +//---------------------------------------------------------------------- + +/** \internal \returns the hyperbolic sine of \a x (coeff-wise). + Uses sinh(x) = (exp(x) - exp(-x)) / 2. + Near overflow, uses sinh(x) = sign(x) * exp(|x|) / 2 via ldexp to avoid inf. + For |x| < 1, uses a direct polynomial to avoid catastrophic cancellation. +*/ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_float(const Packet& x) { + typedef typename unpacket_traits::type Scalar; + static_assert(std::is_same::value, "Scalar type must be float"); + + const Packet sign_mask = pset1(-0.0f); + const Packet abs_x = pandnot(x, sign_mask); + const Packet x_sign = pand(x, sign_mask); + + // For |x| < 1, use a polynomial approximation to avoid + // cancellation in exp(x) - exp(-x). + constexpr float alpha[] = {2.7557314045e-06f, 1.9841270114e-04f, 8.3333335817e-03f, 1.6666666716e-01f}; + const Packet x2 = pmul(x, x); + Packet p_small = ppolevl::run(x2, alpha); + p_small = pmadd(pmul(x2, x), p_small, x); + + // Compute e = exp(|x|) / 2 = exp(|x| - 1) * (e/2), where e is Euler's number. + // Using a single exp avoids a second expensive call, and subtracting 1 (exactly + // representable) instead of ln2 avoids rounding error in the argument to exp, + // which would be amplified into large relative output error. + const Packet half_e = pset1(1.3591409142295225f); // e/2 + const Packet one = pset1(1.0f); + const Packet e = pmul(pexp(psub(abs_x, one)), half_e); + + // Medium path (1 <= |x| < 20): + // sinh(x) = (exp(|x|) - exp(-|x|)) / 2 + // = (2*e - 1/(2*e)) / 2 = e - 1/(4*e) + const Packet quarter = pset1(0.25f); + Packet p_medium = psub(e, pdiv(quarter, e)); + + // Large path (|x| >= 20): exp(-|x|) is negligible, sinh(x) ~ exp(|x|)/2 = e. + const Packet large_threshold = pset1(20.0f); + const Packet large_mask = pcmp_lt(large_threshold, abs_x); + Packet p_large = pselect(large_mask, e, p_medium); + p_large = por(x_sign, p_large); + + const Packet small_mask = pcmp_lt(abs_x, one); + return pselect(small_mask, p_small, p_large); +} + +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_double(const Packet& x) { + typedef typename unpacket_traits::type Scalar; + static_assert(std::is_same::value, "Scalar type must be double"); + + const Packet sign_mask = pset1(-0.0); + const Packet abs_x = pandnot(x, sign_mask); + const Packet x_sign = pand(x, sign_mask); + + // Taylor series: sinh(x) = x + x^3/3! + x^5/5! + ... + x^19/19! + // Polynomial form: sinh(x) = x + x^3 * P(x^2) where P(t) = sum_{k=0}^{8} t^k/(2k+3)! + // ppolevl stores highest-degree coefficient first. + constexpr double alpha[] = { + 8.2206352466243297e-18, // t^8: 1/19! + 2.8114572543455206e-15, // t^7: 1/17! + 7.6471637318198164e-13, // t^6: 1/15! + 1.6059043836821613e-10, // t^5: 1/13! + 2.5052108385441718e-08, // t^4: 1/11! + 2.7557319223985893e-06, // t^3: 1/9! + 1.9841269841269841e-04, // t^2: 1/7! + 8.3333333333333332e-03, // t^1: 1/5! + 1.6666666666666666e-01, // t^0: 1/3! + }; + const Packet x2 = pmul(x, x); + Packet p_small = ppolevl::run(x2, alpha); + p_small = pmadd(pmul(x2, x), p_small, x); + + // Compute e = exp(|x|) / 2 = exp(|x| - 1) * (e/2), where e is Euler's number. + // Subtracting 1 (exactly representable) instead of ln2 avoids rounding error + // in the argument to exp, which would be amplified into large relative error. + const Packet half_e = pset1(1.3591409142295225); // e/2 + const Packet one = pset1(1.0); + const Packet e = pmul(pexp(psub(abs_x, one)), half_e); + + // Medium path (1 <= |x| < 20): + // sinh(x) = (exp(|x|) - exp(-|x|)) / 2 = e - 1/(4*e) + const Packet quarter = pset1(0.25); + Packet p_medium = psub(e, pdiv(quarter, e)); + + // Large path (|x| >= 20): exp(-|x|) is negligible, sinh(x) ~ exp(|x|)/2 = e. + const Packet large_threshold = pset1(20.0); + const Packet large_mask = pcmp_lt(large_threshold, abs_x); + Packet p_large = pselect(large_mask, e, p_medium); + p_large = por(x_sign, p_large); + const Packet small_mask = pcmp_lt(abs_x, one); + return pselect(small_mask, p_small, p_large); +} + +/** \internal \returns the hyperbolic cosine of \a x (coeff-wise). + Uses cosh(x) = (exp(|x|) + exp(-|x|)) / 2. + Near overflow, uses ldexp(exp(|x| - ln2), -1) to avoid premature inf. +*/ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_float(const Packet& x) { + const Packet abs_x = pabs(x); + + // Compute e = exp(|x|) / 2 = exp(|x| - 1) * (e/2), where e is Euler's number. + // Using a single exp avoids a second expensive call, and subtracting 1 (exactly + // representable) instead of ln2 avoids rounding error in the argument to exp, + // which would be amplified into large relative output error. + const Packet half_e = pset1(1.3591409142295225f); // e/2 + const Packet one = pset1(1.0f); + const Packet e = pmul(pexp(psub(abs_x, one)), half_e); + + // Medium path: cosh(x) = (exp(|x|) + exp(-|x|)) / 2 + // = (2*e + 1/(2*e)) / 2 = e + 1/(4*e) + const Packet quarter = pset1(0.25f); + Packet p_medium = padd(e, pdiv(quarter, e)); + + // Large path (|x| >= 20): exp(-|x|) is negligible, cosh(x) ~ exp(|x|)/2 = e. + const Packet large_threshold = pset1(20.0f); + const Packet large_mask = pcmp_lt(large_threshold, abs_x); + return pselect(large_mask, e, p_medium); +} + +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_double(const Packet& x) { + const Packet abs_x = pabs(x); + + // Compute e = exp(|x|) / 2 = exp(|x| - 1) * (e/2), where e is Euler's number. + // Subtracting 1 (exactly representable) instead of ln2 avoids rounding error + // in the argument to exp, which would be amplified into large relative error. + const Packet half_e = pset1(1.3591409142295225); // e/2 + const Packet one = pset1(1.0); + const Packet e = pmul(pexp(psub(abs_x, one)), half_e); + + // Medium path: cosh(x) = (exp(|x|) + exp(-|x|)) / 2 = e + 1/(4*e) + const Packet quarter = pset1(0.25); + Packet p_medium = padd(e, pdiv(quarter, e)); + + // Large path (|x| >= 20): exp(-|x|) is negligible, cosh(x) ~ exp(|x|)/2 = e. + const Packet large_threshold = pset1(20.0); + const Packet large_mask = pcmp_lt(large_threshold, abs_x); + return pselect(large_mask, e, p_medium); +} + +//---------------------------------------------------------------------- +// asinh / acosh +//---------------------------------------------------------------------- + +/** \internal \returns the inverse hyperbolic sine of \a x (coeff-wise). + For small |x|: asinh(x) = sign(x) * log1p(|x| + x^2/(1 + sqrt(1 + x^2))) + For large |x|: asinh(x) = sign(x) * (log(|x|) + ln(2)) + Otherwise: asinh(x) = sign(x) * log(|x| + sqrt(x^2 + 1)) +*/ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_float(const Packet& x) { + const Packet sign_mask = pset1(-0.0f); + const Packet abs_x = pandnot(x, sign_mask); + const Packet x_sign = pand(x, sign_mask); + const Packet one = pset1(1.0f); + + // For |x| < 0.5, use log1p formulation to avoid cancellation: + // asinh(x) = log1p(|x| + x^2 / (1 + sqrt(1 + x^2))) + const Packet x2 = pmul(abs_x, abs_x); + Packet p_small = generic_log1p(padd(abs_x, pdiv(x2, padd(one, psqrt(padd(one, x2)))))); + + // For 0.5 <= |x| < 1e10, use log(|x| + sqrt(x^2 + 1)). + Packet p_med = plog(padd(abs_x, psqrt(padd(x2, one)))); + + // For |x| >= 1e10, use log(2*|x|) = log(|x|) + ln(2) to avoid x^2 overflow. + const Packet ln2 = pset1(0.6931471805599453f); + Packet p_large = padd(plog(abs_x), ln2); + + const Packet small_mask = pcmp_lt(abs_x, pset1(0.5f)); + const Packet large_mask = pcmp_lt(pset1(1e10f), abs_x); + Packet result = pselect(large_mask, p_large, pselect(small_mask, p_small, p_med)); + return por(x_sign, result); +} + +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_double(const Packet& x) { + const Packet sign_mask = pset1(-0.0); + const Packet abs_x = pandnot(x, sign_mask); + const Packet x_sign = pand(x, sign_mask); + const Packet one = pset1(1.0); + + const Packet x2 = pmul(abs_x, abs_x); + Packet p_small = generic_log1p(padd(abs_x, pdiv(x2, padd(one, psqrt(padd(one, x2)))))); + + Packet p_med = plog(padd(abs_x, psqrt(padd(x2, one)))); + + const Packet ln2 = pset1(0.6931471805599453); + Packet p_large = padd(plog(abs_x), ln2); + + const Packet small_mask = pcmp_lt(abs_x, pset1(0.5)); + const Packet large_mask = pcmp_lt(pset1(1e150), abs_x); + Packet result = pselect(large_mask, p_large, pselect(small_mask, p_small, p_med)); + return por(x_sign, result); +} + +/** \internal \returns the inverse hyperbolic cosine of \a x (coeff-wise). + Uses acosh(x) = log(x + sqrt(x^2 - 1)) for x >= 1. + Returns NaN for x < 1. +*/ +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_float(const Packet& x) { + const Packet one = pset1(1.0f); + // For x near 1, use log1p to avoid cancellation: + // acosh(x) = log(x + sqrt(x^2-1)) = log(x + sqrt((x-1)(x+1))) + // For x close to 1, let t = x-1, then: + // acosh(x) = log1p(t + sqrt(t*(t+2))) + const Packet t = psub(x, one); + const Packet small_mask = pcmp_lt(t, pset1(0.5f)); + + // Small path: acosh(x) = log1p(t + sqrt(t*(t+2))) + const Packet two = pset1(2.0f); + Packet p_small = generic_log1p(padd(t, psqrt(pmul(t, padd(t, two))))); + + // Large path: acosh(x) = log(x + sqrt(x^2-1)) + // For very large x, use log(2*x) to avoid overflow in x^2. + const Packet large_threshold = pset1(1e10f); + const Packet huge_mask = pcmp_lt(large_threshold, x); + const Packet x2_safe = pselect(huge_mask, one, pmul(x, x)); + Packet p_large = plog(padd(x, psqrt(psub(x2_safe, one)))); + const Packet log2 = pset1(0.6931471805599453f); + p_large = pselect(huge_mask, padd(plog(x), log2), p_large); + + Packet result = pselect(small_mask, p_small, p_large); + + // Return NaN for x < 1. + const Packet invalid_mask = pcmp_lt(x, one); + return por(invalid_mask, result); +} + +template +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_double(const Packet& x) { + const Packet one = pset1(1.0); + const Packet t = psub(x, one); + const Packet small_mask = pcmp_lt(t, pset1(0.5)); + + // Small path: acosh(x) = log1p(t + sqrt(t*(t+2))) + const Packet two = pset1(2.0); + Packet p_small = generic_log1p(padd(t, psqrt(pmul(t, padd(t, two))))); + + // Large path: acosh(x) = log(x + sqrt(x^2-1)) + const Packet large_threshold = pset1(1e150); + const Packet huge_mask = pcmp_lt(large_threshold, x); + const Packet x2_safe = pselect(huge_mask, one, pmul(x, x)); + Packet p_large = plog(padd(x, psqrt(psub(x2_safe, one)))); + const Packet log2 = pset1(0.6931471805599453); + p_large = pselect(huge_mask, padd(plog(x), log2), p_large); + + Packet result = pselect(small_mask, p_small, p_large); + + const Packet invalid_mask = pcmp_lt(x, one); + return por(invalid_mask, result); +} + } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/arch/Default/Half.h b/Eigen/src/Core/arch/Default/Half.h index 56aeb3fa8..70ee8f900 100644 --- a/Eigen/src/Core/arch/Default/Half.h +++ b/Eigen/src/Core/arch/Default/Half.h @@ -60,12 +60,17 @@ #define EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_F16(PACKET_F, PACKET_F16) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pcos) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, psin) \ + F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, psinh) \ + F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pcosh) \ + F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pasinh) \ + F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pacosh) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pexp) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pexp2) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pexpm1) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog1p) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog2) \ + F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog10) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, preciprocal) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, prsqrt) \ F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pcbrt) \ diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h index 701ddcbee..d7818398a 100644 --- a/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/Eigen/src/Core/arch/NEON/PacketMath.h @@ -202,7 +202,10 @@ struct packet_traits : default_packet_traits { HasASin = 1, HasATan = 1, HasATanh = 1, + HasSinh = 1, + HasCosh = 1, HasLog = 1, + HasLog10 = 1, HasExp = 1, HasLog1p = 1, HasExpm1 = 1, @@ -5040,11 +5043,14 @@ struct packet_traits : default_packet_traits { #if EIGEN_ARCH_ARM64 && !EIGEN_APPLE_DOUBLE_NEON_BUG HasExp = 1, HasLog = 1, + HasLog10 = 1, HasLog1p = 1, HasExpm1 = 1, HasPow = 1, HasATan = 1, HasATanh = 1, + HasSinh = 1, + HasCosh = 1, #endif HasSin = EIGEN_FAST_MATH, HasCos = EIGEN_FAST_MATH, diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 28e635b27..e83240e4f 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -190,8 +190,11 @@ struct packet_traits : default_packet_traits { HasASin = 1, HasATan = 1, HasATanh = 1, + HasSinh = 1, + HasCosh = 1, HasLog = 1, HasLog1p = 1, + HasLog10 = 1, HasExpm1 = 1, HasNdtri = 1, HasExp = 1, @@ -220,10 +223,13 @@ struct packet_traits : default_packet_traits { HasSin = EIGEN_FAST_MATH, HasCos = EIGEN_FAST_MATH, HasTan = EIGEN_FAST_MATH, + HasSinh = 1, + HasCosh = 1, HasTanh = EIGEN_FAST_MATH, HasErf = EIGEN_FAST_MATH, HasErfc = EIGEN_FAST_MATH, HasLog = 1, + HasLog10 = 1, HasExp = 1, HasLog1p = 1, HasExpm1 = 1, diff --git a/benchmarks/Core/bench_cwise_math.cpp b/benchmarks/Core/bench_cwise_math.cpp index 8d7e94e6b..f24b45c2d 100644 --- a/benchmarks/Core/bench_cwise_math.cpp +++ b/benchmarks/Core/bench_cwise_math.cpp @@ -49,8 +49,11 @@ BENCH_CWISE_UNARY(Acos, a.acos(), -0.99, 0.99) BENCH_CWISE_UNARY(Atan, a.atan(), -10, 10) // Hyperbolic / special +BENCH_CWISE_UNARY(Sinh, a.sinh(), -5, 5) +BENCH_CWISE_UNARY(Cosh, a.cosh(), -5, 5) BENCH_CWISE_UNARY(Tanh, a.tanh(), -5, 5) BENCH_CWISE_UNARY(Atanh, a.atanh(), -0.99, 0.99) +BENCH_CWISE_UNARY(Log10, a.log10(), 0.01, 100) BENCH_CWISE_UNARY(Erf, Eigen::erf(a), -4, 4) // Simple operations (should be very fast / memory-bound) @@ -149,8 +152,11 @@ BENCHMARK(BM_Tan) CWISE_SIZES ->Name("Tan_float"); BENCHMARK(BM_Asin) CWISE_SIZES ->Name("Asin_float"); BENCHMARK(BM_Acos) CWISE_SIZES ->Name("Acos_float"); BENCHMARK(BM_Atan) CWISE_SIZES ->Name("Atan_float"); +BENCHMARK(BM_Sinh) CWISE_SIZES ->Name("Sinh_float"); +BENCHMARK(BM_Cosh) CWISE_SIZES ->Name("Cosh_float"); BENCHMARK(BM_Tanh) CWISE_SIZES ->Name("Tanh_float"); BENCHMARK(BM_Atanh) CWISE_SIZES ->Name("Atanh_float"); +BENCHMARK(BM_Log10) CWISE_SIZES ->Name("Log10_float"); BENCHMARK(BM_Erf) CWISE_SIZES ->Name("Erf_float"); BENCHMARK(BM_Abs) CWISE_SIZES ->Name("Abs_float"); BENCHMARK(BM_Square) CWISE_SIZES ->Name("Square_float"); @@ -179,8 +185,11 @@ BENCHMARK(BM_Tan) CWISE_SIZES ->Name("Tan_double"); BENCHMARK(BM_Asin) CWISE_SIZES ->Name("Asin_double"); BENCHMARK(BM_Acos) CWISE_SIZES ->Name("Acos_double"); BENCHMARK(BM_Atan) CWISE_SIZES ->Name("Atan_double"); +BENCHMARK(BM_Sinh) CWISE_SIZES ->Name("Sinh_double"); +BENCHMARK(BM_Cosh) CWISE_SIZES ->Name("Cosh_double"); BENCHMARK(BM_Tanh) CWISE_SIZES ->Name("Tanh_double"); BENCHMARK(BM_Atanh) CWISE_SIZES ->Name("Atanh_double"); +BENCHMARK(BM_Log10) CWISE_SIZES ->Name("Log10_double"); BENCHMARK(BM_Erf) CWISE_SIZES ->Name("Erf_double"); BENCHMARK(BM_Abs) CWISE_SIZES ->Name("Abs_double"); BENCHMARK(BM_Square) CWISE_SIZES ->Name("Square_double");