From e72dfeb8b9fa5662831b5d0bb9d132521f9173dd Mon Sep 17 00:00:00 2001 From: Antonio Sanchez Date: Wed, 3 Mar 2021 09:41:46 -0800 Subject: Fix rint for SSE/NEON. It seems *sometimes* with aggressive optimizations the combination `psub(padd(a, b), b)` trick to force rounding is compiled away. Here we replace with inline assembly to prevent this (I tried `volatile`, but that leads to additional loads from memory). Also fixed an edge case for large inputs `a` where adding `b` bumps the value up a power of two and ends up rounding away more than just the fractional part. If we are over `2^digits` then just return the input. This edge case was missed in the test since the test was comparing approximate equality, which was still satisfied. Adding a strict equality option catches it. --- Eigen/src/Core/arch/NEON/PacketMath.h | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) (limited to 'Eigen/src/Core/arch/NEON/PacketMath.h') diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h index ec6ea90c5..51cebaf2b 100644 --- a/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/Eigen/src/Core/arch/NEON/PacketMath.h @@ -3207,20 +3207,34 @@ template<> EIGEN_STRONG_INLINE Packet4f pceil(const Packet4f& a) template<> EIGEN_STRONG_INLINE Packet4f print(const Packet4f& a) { // Adds and subtracts signum(a) * 2^23 to force rounding. - const Packet4f offset = - pselect(pcmp_lt(a, pzero(a)), - pset1(-static_cast(1<<23)), - pset1(+static_cast(1<<23))); - return psub(padd(a, offset), offset); + const Packet4f limit = pset1(static_cast(1<<23)); + const Packet4f abs_a = pabs(a); + // Inline asm to prevent the compiler from optimizing away the + // addition and subtraction. + // Packet4f r = psub(padd(abs_a, limit), limit); + Packet4f r = abs_a; + __asm__ ("vadd.f32 %[r], %[r], %[limit]\n\t" + "vsub.f32 %[r], %[r], %[limit]" : [r] "+x" (r) : [limit] "x" (limit)); + // If greater than limit, simply return a. Otherwise, account for sign. + r = pselect(pcmp_lt(abs_a, limit), + pselect(pcmp_lt(a, pzero(a)), pnegate(r), r), a); + return r; } template<> EIGEN_STRONG_INLINE Packet2f print(const Packet2f& a) { // Adds and subtracts signum(a) * 2^23 to force rounding. - const Packet2f offset = - pselect(pcmp_lt(a, pzero(a)), - pset1(-static_cast(1<<23)), - pset1(+static_cast(1<<23))); - return psub(padd(a, offset), offset); + const Packet2f limit = pset1(static_cast(1<<23)); + const Packet2f abs_a = pabs(a); + // Inline asm to prevent the compiler from optimizing away the + // addition and subtraction. + // Packet4f r = psub(padd(abs_a, limit), limit); + Packet2f r = abs_a; + __asm__ ("vadd.f32 %[r], %[r], %[limit]\n\t" + "vsub.f32 %[r], %[r], %[limit]" : [r] "+x" (r) : [limit] "x" (limit)); + // If greater than limit, simply return a. Otherwise, account for sign. + r = pselect(pcmp_lt(abs_a, limit), + pselect(pcmp_lt(a, pzero(a)), pnegate(r), r), a); + return r; } template<> EIGEN_STRONG_INLINE Packet4f pfloor(const Packet4f& a) -- cgit v1.2.3