Unify NEON's plog with generic implementation

This commit is contained in:
Gael Guennebaud 2018-11-26 15:02:16 +01:00
parent 2c44c40114
commit c24e98e6a8
2 changed files with 19 additions and 91 deletions

View File

@ -5,13 +5,15 @@
// 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/.
/* The sin, cos, exp, and log functions of this file come from
/* The exp function of this file come from
* Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
*/
#ifndef EIGEN_MATH_FUNCTIONS_NEON_H
#define EIGEN_MATH_FUNCTIONS_NEON_H
#include "../Default/GenericPacketMathFunctions.h"
namespace Eigen {
namespace internal {
@ -87,96 +89,7 @@ Packet4f pexp<Packet4f>(const Packet4f& _x)
template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
Packet4f plog<Packet4f>(const Packet4f& _x)
{
Packet4f x = _x;
_EIGEN_DECLARE_CONST_Packet4f(1 , 1.0f);
_EIGEN_DECLARE_CONST_Packet4f(half, 0.5f);
_EIGEN_DECLARE_CONST_Packet4i(0x7f, 0x7f);
const Packet4f p4f_minus_inf = vreinterpretq_f32_s32(pset1<Packet4i>(0xff800000));
_EIGEN_DECLARE_CONST_Packet4i(inv_mant_mask, ~0x7f800000);
/* natural logarithm computed for 4 simultaneous float
return NaN for x <= 0
*/
_EIGEN_DECLARE_CONST_Packet4f(cephes_SQRTHF, 0.707106781186547524f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p0, 7.0376836292E-2f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p1, - 1.1514610310E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p2, 1.1676998740E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p3, - 1.2420140846E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p4, + 1.4249322787E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p5, - 1.6668057665E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p6, + 2.0000714765E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p7, - 2.4999993993E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_p8, + 3.3333331174E-1f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_q1, -2.12194440e-4f);
_EIGEN_DECLARE_CONST_Packet4f(cephes_log_q2, 0.693359375f);
x = vmaxq_f32(x, vdupq_n_f32(0)); /* force flush to zero on denormal values */
Packet4f iszero_mask = vreinterpretq_f32_u32(vceqq_f32(_x, vdupq_n_f32(0)));
Packet4f invalid_mask = vreinterpretq_f32_u32(vmvnq_u32(vcgeq_f32(_x, vdupq_n_f32(0))));
Packet4i ux = vreinterpretq_s32_f32(x);
Packet4i emm0 = vshrq_n_s32(ux, 23);
/* keep only the fractional part */
ux = vandq_s32(ux, p4i_inv_mant_mask);
ux = vorrq_s32(ux, vreinterpretq_s32_f32(p4f_half));
x = vreinterpretq_f32_s32(ux);
emm0 = vsubq_s32(emm0, p4i_0x7f);
Packet4f e = vcvtq_f32_s32(emm0);
e = vaddq_f32(e, p4f_1);
/* part2:
if( x < SQRTHF ) {
e -= 1;
x = x + x - 1.0;
} else { x = x - 1.0; }
*/
Packet4ui mask = vcltq_f32(x, p4f_cephes_SQRTHF);
Packet4f tmp = vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(x), mask));
x = vsubq_f32(x, p4f_1);
e = vsubq_f32(e, vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(p4f_1), mask)));
x = vaddq_f32(x, tmp);
Packet4f z = vmulq_f32(x,x);
Packet4f y = p4f_cephes_log_p0;
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p1);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p2);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p3);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p4);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p5);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p6);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p7);
y = vmulq_f32(y, x);
y = vaddq_f32(y, p4f_cephes_log_p8);
y = vmulq_f32(y, x);
y = vmulq_f32(y, z);
tmp = vmulq_f32(e, p4f_cephes_log_q1);
y = vaddq_f32(y, tmp);
tmp = vmulq_f32(z, p4f_half);
y = vsubq_f32(y, tmp);
tmp = vmulq_f32(e, p4f_cephes_log_q2);
x = vaddq_f32(x, y);
x = vaddq_f32(x, tmp);
x = por(x, invalid_mask);
x = por(pandnot(x,iszero_mask), pand(iszero_mask, p4f_minus_inf));
return x;
return plog_float(_x);
}
} // end namespace internal

View File

@ -145,6 +145,8 @@ template<> struct unpacket_traits<Packet4i> { typedef int32_t type; enum {size=4
template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float& from) { return vdupq_n_f32(from); }
template<> EIGEN_STRONG_INLINE Packet4i pset1<Packet4i>(const int32_t& from) { return vdupq_n_s32(from); }
template<> EIGEN_STRONG_INLINE Packet4f pset1frombits<Packet4f>(unsigned int from) { return vreinterpretq_f32_u32(vdupq_n_u32(from)); }
template<> EIGEN_STRONG_INLINE Packet4f plset<Packet4f>(const float& a)
{
const float f[] = {0, 1, 2, 3};
@ -249,6 +251,11 @@ template<> EIGEN_STRONG_INLINE Packet4i pmin<Packet4i>(const Packet4i& a, const
template<> EIGEN_STRONG_INLINE Packet4f pmax<Packet4f>(const Packet4f& a, const Packet4f& b) { return vmaxq_f32(a,b); }
template<> EIGEN_STRONG_INLINE Packet4i pmax<Packet4i>(const Packet4i& a, const Packet4i& b) { return vmaxq_s32(a,b); }
template<> EIGEN_STRONG_INLINE Packet4f pcmp_le(const Packet4f& a, const Packet4f& b) { return vreinterpretq_f32_u32(vcleq_f32(a,b)); }
template<> EIGEN_STRONG_INLINE Packet4f pcmp_lt(const Packet4f& a, const Packet4f& b) { return vreinterpretq_f32_u32(vcltq_f32(a,b)); }
template<> EIGEN_STRONG_INLINE Packet4f pcmp_eq(const Packet4f& a, const Packet4f& b) { return vreinterpretq_f32_u32(vceqq_f32(a,b)); }
template<> EIGEN_STRONG_INLINE Packet4f pcmp_lt_or_nan(const Packet4f& a, const Packet4f& b) { return vreinterpretq_f32_u32(vmvnq_u32(vcgeq_f32(a,b))); }
// Logical Operations are not supported for float, so we have to reinterpret casts using NEON intrinsics
template<> EIGEN_STRONG_INLINE Packet4f pand<Packet4f>(const Packet4f& a, const Packet4f& b)
{
@ -364,6 +371,14 @@ template<> EIGEN_STRONG_INLINE Packet4i preverse(const Packet4i& a) {
template<> EIGEN_STRONG_INLINE Packet4f pabs(const Packet4f& a) { return vabsq_f32(a); }
template<> EIGEN_STRONG_INLINE Packet4i pabs(const Packet4i& a) { return vabsq_s32(a); }
template<> EIGEN_STRONG_INLINE Packet4f pshiftright_and_cast(Packet4f a, int n) {
return vcvtq_f32_s32(vshrq_n_s32(vreinterpretq_s32_f32(a),n));
}
template<> EIGEN_STRONG_INLINE Packet4f pfrexp<Packet4f>(const Packet4f& a, Packet4f& exponent) {
return pfrexp_float(a,exponent);
}
template<> EIGEN_STRONG_INLINE float predux<Packet4f>(const Packet4f& a)
{
float32x2_t a_lo, a_hi, sum;