Unify NEON's plog with generic implementation

2 files changed, 19 insertions, 91 deletions

diff --git a/Eigen/src/Core/arch/NEON/MathFunctions.h b/Eigen/src/Core/arch/NEON/MathFunctions.h
index d218c8851..8a70b771c 100644
--- a/Eigen/src/Core/arch/NEON/MathFunctions.h
+++ b/Eigen/src/Core/arch/NEON/MathFunctions.h
@@ -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
diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h
index 010739380..72f076e50 100644
--- a/Eigen/src/Core/arch/NEON/PacketMath.h
+++ b/Eigen/src/Core/arch/NEON/PacketMath.h
@@ -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;