devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Vectorize sinh, cosh, and log10

Browse Source 
libeigen/eigen!2368

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-04-01 20:41:18 -07:00
parent 30e669cfe1
commit 9513d3878e
12 changed files with 393 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
Eigen/src/Core/GenericPacketMath.h
View File

@@ -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 <typename Packet>
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 <typename Packet>
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 <typename Packet>
EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp(const Packet& a) {

5
Eigen/src/Core/arch/AVX/MathFunctions.h
View File

@@ -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)

6
Eigen/src/Core/arch/AVX/PacketMath.h
View File

@@ -117,7 +117,10 @@ struct packet_traits<float> : 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<double> : 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,

6
Eigen/src/Core/arch/AVX512/PacketMath.h
View File

@@ -172,10 +172,13 @@ struct packet_traits<float> : 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<double> : 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,

5
Eigen/src/Core/arch/Default/BFloat16.h
View File

@@ -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) \

19
Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
View File

@@ -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 <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_float(const Packet& x) {
const Packet cst_log10e = pset1<Packet>(0.4342944819032518f);
return pmul(plog(x), cst_log10e);
}
/** \internal \returns log10(x) for double precision float.
Computed as log(x) * log10(e). */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_double(const Packet& x) {
const Packet cst_log10e = pset1<Packet>(0.4342944819032518);
return pmul(plog(x), cst_log10e);
}
} // end namespace internal
} // end namespace Eigen

54
Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h
View File

@@ -150,6 +150,46 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_float(const Pa
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_double(const Packet& x);
/** \internal \returns sinh(x) for single precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_float(const Packet& x);
/** \internal \returns sinh(x) for double precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_double(const Packet& x);
/** \internal \returns cosh(x) for single precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_float(const Packet& x);
/** \internal \returns cosh(x) for double precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_double(const Packet& x);
/** \internal \returns asinh(x) for single precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_float(const Packet& x);
/** \internal \returns asinh(x) for double precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_double(const Packet& x);
/** \internal \returns acosh(x) for single precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_float(const Packet& x);
/** \internal \returns acosh(x) for double precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_double(const Packet& x);
/** \internal \returns log10(x) for single precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_float(const Packet& x);
/** \internal \returns log10(x) for double precision float */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_double(const Packet& x);
/** \internal \returns sqrt(x) for complex types */
template <typename Packet>
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) \

258
Eigen/src/Core/arch/Default/GenericPacketMathTrig.h
View File

@@ -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 <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_float(const Packet& x) {
typedef typename unpacket_traits<Packet>::type Scalar;
static_assert(std::is_same<Scalar, float>::value, "Scalar type must be float");
const Packet sign_mask = pset1<Packet>(-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<Packet, 3>::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<Packet>(1.3591409142295225f); // e/2
const Packet one = pset1<Packet>(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<Packet>(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<Packet>(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 <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_double(const Packet& x) {
typedef typename unpacket_traits<Packet>::type Scalar;
static_assert(std::is_same<Scalar, double>::value, "Scalar type must be double");
const Packet sign_mask = pset1<Packet>(-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<Packet, 8>::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<Packet>(1.3591409142295225); // e/2
const Packet one = pset1<Packet>(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<Packet>(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<Packet>(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 <typename Packet>
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<Packet>(1.3591409142295225f); // e/2
const Packet one = pset1<Packet>(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<Packet>(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<Packet>(20.0f);
const Packet large_mask = pcmp_lt(large_threshold, abs_x);
return pselect(large_mask, e, p_medium);
}
template <typename Packet>
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<Packet>(1.3591409142295225); // e/2
const Packet one = pset1<Packet>(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<Packet>(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<Packet>(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 <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_float(const Packet& x) {
const Packet sign_mask = pset1<Packet>(-0.0f);
const Packet abs_x = pandnot(x, sign_mask);
const Packet x_sign = pand(x, sign_mask);
const Packet one = pset1<Packet>(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<Packet>(0.6931471805599453f);
Packet p_large = padd(plog(abs_x), ln2);
const Packet small_mask = pcmp_lt(abs_x, pset1<Packet>(0.5f));
const Packet large_mask = pcmp_lt(pset1<Packet>(1e10f), abs_x);
Packet result = pselect(large_mask, p_large, pselect(small_mask, p_small, p_med));
return por(x_sign, result);
}
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_double(const Packet& x) {
const Packet sign_mask = pset1<Packet>(-0.0);
const Packet abs_x = pandnot(x, sign_mask);
const Packet x_sign = pand(x, sign_mask);
const Packet one = pset1<Packet>(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<Packet>(0.6931471805599453);
Packet p_large = padd(plog(abs_x), ln2);
const Packet small_mask = pcmp_lt(abs_x, pset1<Packet>(0.5));
const Packet large_mask = pcmp_lt(pset1<Packet>(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 <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_float(const Packet& x) {
const Packet one = pset1<Packet>(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<Packet>(0.5f));
// Small path: acosh(x) = log1p(t + sqrt(t*(t+2)))
const Packet two = pset1<Packet>(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<Packet>(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<Packet>(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 <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pacosh_double(const Packet& x) {
const Packet one = pset1<Packet>(1.0);
const Packet t = psub(x, one);
const Packet small_mask = pcmp_lt(t, pset1<Packet>(0.5));
// Small path: acosh(x) = log1p(t + sqrt(t*(t+2)))
const Packet two = pset1<Packet>(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<Packet>(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<Packet>(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

5
Eigen/src/Core/arch/Default/Half.h
View File

@@ -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) \

6
Eigen/src/Core/arch/NEON/PacketMath.h
View File

@@ -202,7 +202,10 @@ struct packet_traits<float> : 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<double> : 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,

6
Eigen/src/Core/arch/SSE/PacketMath.h
View File

@@ -190,8 +190,11 @@ struct packet_traits<float> : 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<double> : 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,