Core bug fixes: MSA, AltiVec, RVV, LSX and SVE backends libeigen/eigen!3349 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index 3a645a1..f21dfb6 100644 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h
@@ -1136,6 +1136,8 @@ // Do one Newton-Raphson iteration to get the needed accuracy t = vec_nmsub(y_0, b, p4f_ONE); y_1 = vec_madd(y_0, t, y_0); + // For b = +-0 or +-inf the step computes 0 * inf = NaN, while the estimate +-inf or +-0 is exact. + y_1 = vec_sel(y_0, y_1, vec_cmpeq(t, t)); return vec_madd(a, y_1, p4f_MZERO); #else @@ -1261,8 +1263,9 @@ EIGEN_STRONG_INLINE Packet4f pcmp_le(const Packet4f& a, const Packet4f& b) { return reinterpret_cast<Packet4f>(vec_cmple(a, b)); } -// To fix bug with vec_cmplt on older versions -#ifdef EIGEN_VECTORIZE_VSX +// To fix bug with vec_cmplt on older versions. Clang needs it: for AltiVec vectors its a < b is one bool, true when +// every lane compares less, which breaks the generic pcmp_lt. +#if defined(EIGEN_VECTORIZE_VSX) || EIGEN_COMP_CLANG template <> EIGEN_STRONG_INLINE Packet4f pcmp_lt(const Packet4f& a, const Packet4f& b) { return reinterpret_cast<Packet4f>(vec_cmplt(a, b));
diff --git a/Eigen/src/Core/arch/MSA/Complex.h b/Eigen/src/Core/arch/MSA/Complex.h index d35884f..bf9a637 100644 --- a/Eigen/src/Core/arch/MSA/Complex.h +++ b/Eigen/src/Core/arch/MSA/Complex.h
@@ -103,9 +103,6 @@ HasMul = 1, HasDiv = 1, HasNegate = 1, - HasSqrt = 1, - HasLog = 1, - HasExp = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -400,8 +397,6 @@ HasMul = 1, HasDiv = 1, HasNegate = 1, - HasSqrt = 1, - HasLog = 1, HasAbs = 0, HasAbs2 = 0, HasMin = 0, @@ -595,9 +590,6 @@ kernel.packet[1].v = v2; } -EIGEN_INSTANTIATE_COMPLEX_MATH_FUNCS(Packet2cf) -EIGEN_INSTANTIATE_COMPLEX_MATH_FUNCS_NO_EXP(Packet1cd) - EIGEN_GCC_FAST_MATH_COMPLEX_VECTORIZE_WORKAROUND_POP } // end namespace internal
diff --git a/Eigen/src/Core/arch/MSA/MathFunctions.h b/Eigen/src/Core/arch/MSA/MathFunctions.h index 947453f..fef24ec 100644 --- a/Eigen/src/Core/arch/MSA/MathFunctions.h +++ b/Eigen/src/Core/arch/MSA/MathFunctions.h
@@ -15,7 +15,7 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. // SPDX-License-Identifier: MPL-2.0 -/* The sin, cos, exp, and log functions of this file come from +/* The exp and log functions of this file come from * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/ */ @@ -225,96 +225,18 @@ return p; } -template <bool sine> -Packet4f psincos_inner_msa_float(const Packet4f& _x) { - static EIGEN_DECLARE_CONST_Packet4f(sincos_max_arg, 13176795.0f); // Approx. (2**24) / (4/Pi). - static EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP1, -0.78515625f); - static EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP2, -2.4187564849853515625e-4f); - static EIGEN_DECLARE_CONST_Packet4f(minus_cephes_DP3, -3.77489497744594108e-8f); - static EIGEN_DECLARE_CONST_Packet4f(sincof_p0, -1.9515295891e-4f); - static EIGEN_DECLARE_CONST_Packet4f(sincof_p1, 8.3321608736e-3f); - static EIGEN_DECLARE_CONST_Packet4f(sincof_p2, -1.6666654611e-1f); - static EIGEN_DECLARE_CONST_Packet4f(coscof_p0, 2.443315711809948e-5f); - static EIGEN_DECLARE_CONST_Packet4f(coscof_p1, -1.388731625493765e-3f); - static EIGEN_DECLARE_CONST_Packet4f(coscof_p2, 4.166664568298827e-2f); - static EIGEN_DECLARE_CONST_Packet4f(cephes_FOPI, 1.27323954473516f); // 4/Pi. - static EIGEN_DECLARE_CONST_Packet4f(half, 0.5f); - static EIGEN_DECLARE_CONST_Packet4f(1, 1.0f); +EIGEN_FLOAT_PACKET_FUNCTION(sin, Packet4f) +EIGEN_FLOAT_PACKET_FUNCTION(cos, Packet4f) - Packet4f x = pabs(_x); - - // Translate infinite arguments into NANs. - Packet4f zero_or_nan_if_inf = psub(_x, _x); - x = padd(x, zero_or_nan_if_inf); - // Prevent sin/cos from generating values larger than 1.0 in magnitude - // for very large arguments by setting x to 0.0. - Packet4i small_or_nan_mask = __builtin_msa_fcult_w(x, p4f_sincos_max_arg); - x = pand(x, (Packet4f)small_or_nan_mask); - - // Scale x by 4/Pi to find x's octant. - Packet4f y = pmul(x, p4f_cephes_FOPI); - // Get the octant. We'll reduce x by this number of octants or by one more than it. - Packet4i y_int = __builtin_msa_ftrunc_s_w(y); - // x's from even-numbered octants will translate to octant 0: [0, +Pi/4]. - // x's from odd-numbered octants will translate to octant -1: [-Pi/4, 0]. - // Adjustment for odd-numbered octants: octant = (octant + 1) & (~1). - Packet4i y_int1 = __builtin_msa_addvi_w(y_int, 1); - Packet4i y_int2 = (Packet4i)__builtin_msa_bclri_w((Packet4ui)y_int1, 0); // bclri = bit-clear - y = __builtin_msa_ffint_s_w(y_int2); - - // Compute the sign to apply to the polynomial. - Packet4i sign_mask = sine ? pxor(__builtin_msa_slli_w(y_int1, 29), (Packet4i)_x) - : __builtin_msa_slli_w(__builtin_msa_addvi_w(y_int, 3), 29); - - // Get the polynomial selection mask. - // We'll calculate both (sin and cos) polynomials and then select from the two. - Packet4i poly_mask = __builtin_msa_ceqi_w(__builtin_msa_slli_w(y_int2, 30), 0); - - // Reduce x by y octants to get: -Pi/4 <= x <= +Pi/4. - // The magic pass: "Extended precision modular arithmetic" - // x = ((x - y * DP1) - y * DP2) - y * DP3 - Packet4f tmp1 = pmul(y, p4f_minus_cephes_DP1); - Packet4f tmp2 = pmul(y, p4f_minus_cephes_DP2); - Packet4f tmp3 = pmul(y, p4f_minus_cephes_DP3); - x = padd(x, tmp1); - x = padd(x, tmp2); - x = padd(x, tmp3); - - // Evaluate the cos(x) polynomial. - y = p4f_coscof_p0; - Packet4f z = pmul(x, x); - y = pmadd(y, z, p4f_coscof_p1); - y = pmadd(y, z, p4f_coscof_p2); - y = pmul(y, z); - y = pmul(y, z); - y = __builtin_msa_fmsub_w(y, z, p4f_half); - y = padd(y, p4f_1); - - // Evaluate the sin(x) polynomial. - Packet4f y2 = p4f_sincof_p0; - y2 = pmadd(y2, z, p4f_sincof_p1); - y2 = pmadd(y2, z, p4f_sincof_p2); - y2 = pmul(y2, z); - y2 = pmadd(y2, x, x); - - // Select the correct result from the two polynomials. - y = sine ? (Packet4f)__builtin_msa_bsel_v((v16u8)poly_mask, (v16u8)y, (v16u8)y2) - : (Packet4f)__builtin_msa_bsel_v((v16u8)poly_mask, (v16u8)y2, (v16u8)y); - - // Update the sign. - sign_mask = pxor(sign_mask, (Packet4i)y); - y = (Packet4f)__builtin_msa_binsli_w((v4u32)y, (v4u32)sign_mask, 0); // binsli = bit-insert-left - return y; +// generic_exp2 clamps its argument with pmin and pmax, which keep the number over a NaN under EIGEN_FAST_MATH. +template <> +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet4f pexp2<Packet4f>(const Packet4f& x) { + return pselect(pisnan(x), x, generic_exp2(x)); } template <> -EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet4f psin<Packet4f>(const Packet4f& x) { - return psincos_inner_msa_float</* sine */ true>(x); -} - -template <> -EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet4f pcos<Packet4f>(const Packet4f& x) { - return psincos_inner_msa_float</* sine */ false>(x); +EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet2d pexp2<Packet2d>(const Packet2d& x) { + return pselect(pisnan(x), x, generic_exp2(x)); } template <>
diff --git a/Eigen/src/Core/arch/MSA/PacketMath.h b/Eigen/src/Core/arch/MSA/PacketMath.h index cce44b9..0a689b8 100644 --- a/Eigen/src/Core/arch/MSA/PacketMath.h +++ b/Eigen/src/Core/arch/MSA/PacketMath.h
@@ -108,6 +108,7 @@ masked_store_available = false }; typedef Packet4f half; + typedef Packet4i integer_packet; }; template <> @@ -544,6 +545,31 @@ return __builtin_msa_add_a_w(zero, a); } +template <int N> +EIGEN_STRONG_INLINE Packet4i parithmetic_shift_right(const Packet4i& a) { + return __builtin_msa_srai_w(a, N); +} + +template <int N> +EIGEN_STRONG_INLINE Packet4i plogical_shift_right(const Packet4i& a) { + return __builtin_msa_srli_w(a, N); +} + +template <int N> +EIGEN_STRONG_INLINE Packet4i plogical_shift_left(const Packet4i& a) { + return __builtin_msa_slli_w(a, N); +} + +template <> +EIGEN_STRONG_INLINE Packet4i pcast<Packet4f, Packet4i>(const Packet4f& a) { + return __builtin_msa_ftrunc_s_w(a); +} + +template <> +EIGEN_STRONG_INLINE Packet4f pldexp<Packet4f>(const Packet4f& a, const Packet4f& exponent) { + return pldexp_generic(a, exponent); +} + template <> EIGEN_STRONG_INLINE float predux<Packet4f>(const Packet4f& a) { EIGEN_MSA_DEBUG; @@ -767,26 +793,6 @@ } template <> -EIGEN_STRONG_INLINE Packet4f pround<Packet4f>(const Packet4f& a) { - Packet4f v = a; - int32_t old_mode, new_mode; - asm volatile( - "cfcmsa %[old_mode], $1\n" - "ori %[new_mode], %[old_mode], 3\n" - "xori %[new_mode], %[new_mode], 3\n" // 0 = round to nearest, ties to even. - "ctcmsa $1, %[new_mode]\n" - "frint.w %w[v], %w[v]\n" - "ctcmsa $1, %[old_mode]\n" - : // outputs - [old_mode] "=r"(old_mode), [new_mode] "=r"(new_mode), - [v] "+f"(v) - : // inputs - : // clobbers - ); - return v; -} - -template <> EIGEN_STRONG_INLINE Packet4f print<Packet4f>(const Packet4f& a) { // frint.w uses the current rounding mode (default: round to nearest, ties to even). Packet4f v = a; @@ -871,9 +877,63 @@ masked_store_available = false }; typedef Packet2d half; + typedef Packet2l integer_packet; }; template <> +struct unpacket_traits<Packet2l> { + typedef int64_t type; + enum { + size = 2, + alignment = Aligned16, + vectorizable = true, + masked_load_available = false, + masked_store_available = false + }; + typedef Packet2l half; +}; + +template <> +EIGEN_STRONG_INLINE Packet2l pset1<Packet2l>(const int64_t& from) { + return __builtin_msa_fill_d(from); +} + +template <int N> +EIGEN_STRONG_INLINE Packet2l plogical_shift_right(const Packet2l& a) { + return __builtin_msa_srli_d(a, N); +} + +template <int N> +EIGEN_STRONG_INLINE Packet2l plogical_shift_left(const Packet2l& a) { + return __builtin_msa_slli_d(a, N); +} + +template <> +EIGEN_STRONG_INLINE Packet2l pcast<Packet2d, Packet2l>(const Packet2d& a) { + return __builtin_msa_ftrunc_s_d(a); +} + +template <> +EIGEN_STRONG_INLINE Packet2l pand<Packet2l>(const Packet2l& a, const Packet2l& b) { + return (Packet2l)__builtin_msa_and_v((v16u8)a, (v16u8)b); +} + +template <> +EIGEN_STRONG_INLINE Packet2l por<Packet2l>(const Packet2l& a, const Packet2l& b) { + return (Packet2l)__builtin_msa_or_v((v16u8)a, (v16u8)b); +} + +template <> +EIGEN_STRONG_INLINE Packet2l pxor<Packet2l>(const Packet2l& a, const Packet2l& b) { + return (Packet2l)__builtin_msa_xor_v((v16u8)a, (v16u8)b); +} + +template <> +EIGEN_STRONG_INLINE Packet2l pandnot<Packet2l>(const Packet2l& a, const Packet2l& b) { + return pand(a, (Packet2l)__builtin_msa_xori_b((v16u8)b, 255)); +} + +template <> EIGEN_STRONG_INLINE Packet2d pset1<Packet2d>(const double& from) { EIGEN_MSA_DEBUG; @@ -1201,26 +1261,6 @@ } template <> -EIGEN_STRONG_INLINE Packet2d pround<Packet2d>(const Packet2d& a) { - Packet2d v = a; - int32_t old_mode, new_mode; - asm volatile( - "cfcmsa %[old_mode], $1\n" - "ori %[new_mode], %[old_mode], 3\n" - "xori %[new_mode], %[new_mode], 3\n" // 0 = round to nearest, ties to even. - "ctcmsa $1, %[new_mode]\n" - "frint.d %w[v], %w[v]\n" - "ctcmsa $1, %[old_mode]\n" - : // outputs - [old_mode] "=r"(old_mode), [new_mode] "=r"(new_mode), - [v] "+f"(v) - : // inputs - : // clobbers - ); - return v; -} - -template <> EIGEN_STRONG_INLINE Packet2d print<Packet2d>(const Packet2d& a) { // frint.d uses the current rounding mode (default: round to nearest, ties to even). Packet2d v = a; @@ -1253,6 +1293,11 @@ return (Packet2d)__builtin_msa_fcult_d(a, b); } +template <> +EIGEN_STRONG_INLINE Packet2d pldexp<Packet2d>(const Packet2d& a, const Packet2d& exponent) { + return pldexp_generic(a, exponent); +} + } // end namespace internal } // end namespace Eigen
diff --git a/Eigen/src/Core/arch/SVE/PacketMath.h b/Eigen/src/Core/arch/SVE/PacketMath.h index 87381c9..a6120dd 100644 --- a/Eigen/src/Core/arch/SVE/PacketMath.h +++ b/Eigen/src/Core/arch/SVE/PacketMath.h
@@ -261,18 +261,29 @@ EIGEN_DEBUG_UNALIGNED_STORE svst1_s32(svptrue_b32(), to, from); } +// The 32-bit gathers and scatters address lane i through the int32 index i * stride. Where that would overflow, the +// lanes are addressed one at a time. +EIGEN_STRONG_INLINE bool sve_s32index_fits(Index stride) { + const Index max_stride = Index(NumTraits<numext::int32_t>::highest()) / Index(EIGEN_ARM64_SVE_VL / 32 - 1); + return stride <= max_stride && stride >= -max_stride; +} + template <> EIGEN_DEVICE_FUNC inline PacketXi pgather<numext::int32_t, PacketXi>(const numext::int32_t* from, Index stride) { + if (!sve_s32index_fits(stride)) + return pgather_partial<numext::int32_t, PacketXi>(from, stride, unpacket_traits<PacketXi>::size); // Index format: {base=0, base+stride, base+stride*2, base+stride*3, ...} - svint32_t indices = svindex_s32(0, stride); + svint32_t indices = svindex_s32(0, numext::int32_t(stride)); return svld1_gather_s32index_s32(svptrue_b32(), from, indices); } template <> EIGEN_DEVICE_FUNC inline void pscatter<numext::int32_t, PacketXi>(numext::int32_t* to, const PacketXi& from, Index stride) { + if (!sve_s32index_fits(stride)) + return pscatter_partial<numext::int32_t, PacketXi>(to, from, stride, unpacket_traits<PacketXi>::size); // Index format: {base=0, base+stride, base+stride*2, base+stride*3, ...} - svint32_t indices = svindex_s32(0, stride); + svint32_t indices = svindex_s32(0, numext::int32_t(stride)); svst1_scatter_s32index_s32(svptrue_b32(), to, indices, from); } @@ -888,15 +899,19 @@ template <> EIGEN_DEVICE_FUNC inline PacketXf pgather<float, PacketXf>(const float* from, Index stride) { + if (!sve_s32index_fits(stride)) + return pgather_partial<float, PacketXf>(from, stride, unpacket_traits<PacketXf>::size); // Index format: {base=0, base+stride, base+stride*2, base+stride*3, ...} - svint32_t indices = svindex_s32(0, stride); + svint32_t indices = svindex_s32(0, numext::int32_t(stride)); return svld1_gather_s32index_f32(svptrue_b32(), from, indices); } template <> EIGEN_DEVICE_FUNC inline void pscatter<float, PacketXf>(float* to, const PacketXf& from, Index stride) { + if (!sve_s32index_fits(stride)) + return pscatter_partial<float, PacketXf>(to, from, stride, unpacket_traits<PacketXf>::size); // Index format: {base=0, base+stride, base+stride*2, base+stride*3, ...} - svint32_t indices = svindex_s32(0, stride); + svint32_t indices = svindex_s32(0, numext::int32_t(stride)); svst1_scatter_s32index_f32(svptrue_b32(), to, indices, from); }
diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp index 18f5069..a12d541 100644 --- a/test/array_cwise.cpp +++ b/test/array_cwise.cpp
@@ -960,6 +960,54 @@ VERIFY_IS_APPROX(m1.cosh().acosh(), acosh(cosh(m1))); VERIFY_IS_APPROX(m1.logistic(), logistic(m1)); + // sin and cos must reduce large arguments accurately. A Cody-Waite reduction in Scalar loses every bit of the + // reduced argument long before |x| reaches 2^digits. + if (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value) { + Array<Scalar, Dynamic, 1> x(2 * NumTraits<Scalar>::digits()); + for (Index i = 0; i < x.size(); ++i) x(i) = Scalar(std::ldexp(1.0 + 0.37 * double(i % 3), int(i))); + const Array<Scalar, Dynamic, 1> s = x.sin(), c = x.cos(); + for (Index i = 0; i < x.size(); ++i) { + VERIFY(abs(s(i) - numext::sin(x(i))) <= Scalar(4) * NumTraits<Scalar>::epsilon()); + VERIFY(abs(c(i) - numext::cos(x(i))) <= Scalar(4) * NumTraits<Scalar>::epsilon()); + } + // tanh(x) rounds to +-1 for |x| >= 20. The rational approximation clamps |x| where it reaches 1, which depends on + // whether pmadd is fused. + const Scalar inf = NumTraits<Scalar>::infinity(); + const Scalar big[] = {inf, -inf, Scalar(20), Scalar(-20)}; + for (Index i = 0; i < x.size(); ++i) x(i) = big[i % 4]; + // Evaluate into an array: comparing the expression would call the scalar tanh, which is std::tanh for double. + const Array<Scalar, Dynamic, 1> t = x.tanh(), sgn = x.sign(); + VERIFY_IS_CWISE_EQUAL(t, sgn); + } + + // Division by zero and by infinity is exact, also where pdiv refines a reciprocal estimate. + { + const Scalar inf = NumTraits<Scalar>::infinity(); + Array<Scalar, Dynamic, 1> num(32), den(32), quotient(32); + for (Index i = 0; i < num.size(); ++i) { + num(i) = (i % 4 < 2) ? Scalar(3) : Scalar(-3); + den(i) = (i % 2) ? inf : Scalar(0); + quotient(i) = (i % 2) ? Scalar(0) : num(i) * inf; + } + // Evaluate into an array: comparing the expression would divide coefficient by coefficient. + const Array<Scalar, Dynamic, 1> q = num / den; + VERIFY_IS_CWISE_EQUAL(q, quotient); + } + + // rint rounds halfway cases to even. generic_rint passes values beyond 2^digits through by a lane-wise pcmp_lt, so + // every packet mixes both kinds. + { + const float halves[] = {0.5f, 1.5f, 2.5f, 1e30f, -0.5f, -1.5f, -2.5f, -1e30f}; + const float even[] = {0.f, 2.f, 2.f, 1e30f, -0.f, -2.f, -2.f, -1e30f}; + Array<Scalar, Dynamic, 1> h(32), r(32); + for (Index i = 0; i < h.size(); ++i) { + h(i) = Scalar(halves[i % 8]); + r(i) = Scalar(even[i % 8]); + } + const Array<Scalar, Dynamic, 1> rh = h.rint(); + VERIFY_IS_CWISE_EQUAL(rh, r); + } + VERIFY_IS_APPROX(m1.arg(), arg(m1)); // carg of a real array is its arg as a real Scalar, on the vectorized path too. VERIFY_IS_CWISE_EQUAL(m1.carg(), m1.arg());