Merged in rmlarsen/eigen (pull request PR-230)

Fix a bug in psqrt for SSE and AVX when EIGEN_FAST_MATH=1

3 files changed, 48 insertions, 46 deletions

diff --git a/Eigen/src/Core/arch/AVX/MathFunctions.h b/Eigen/src/Core/arch/AVX/MathFunctions.h
index d21ec39dd..6af67ce2d 100644
--- a/Eigen/src/Core/arch/AVX/MathFunctions.h
+++ b/Eigen/src/Core/arch/AVX/MathFunctions.h
@@ -355,30 +355,27 @@ pexp<Packet4d>(const Packet4d& _x) {
 // Functions for sqrt.
 // The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step
 // of Newton's method, at a cost of 1-2 bits of precision as opposed to the
-// exact solution. The main advantage of this approach is not just speed, but
-// also the fact that it can be inlined and pipelined with other computations,
-// further reducing its effective latency.
+// exact solution. It does not handle +inf, or denormalized numbers correctly.
+// The main advantage of this approach is not just speed, but also the fact that
+// it can be inlined and pipelined with other computations, further reducing its
+// effective latency. This is similar to Quake3's fast inverse square root.
+// For detail see here: http://www.beyond3d.com/content/articles/8/
 #if EIGEN_FAST_MATH
 template <>
 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f
 psqrt<Packet8f>(const Packet8f& _x) {
-  _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f);
-  _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f);
-  _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000);
-
-  Packet8f neg_half = pmul(_x, p8f_minus_half);
-
-  // select only the inverse sqrt of positive normal inputs (denormals are
-  // flushed to zero and cause infs as well).
-  Packet8f non_zero_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_GE_OQ);
-  Packet8f x = _mm256_and_ps(non_zero_mask, _mm256_rsqrt_ps(_x));
-
+  Packet8f half = pmul(_x, pset1<Packet8f>(.5f));
+  Packet8f denormal_mask = _mm256_and_ps(
+      _mm256_cmp_ps(_x, pset1<Packet8f>((std::numeric_limits<float>::min)()),
+                    _CMP_LT_OQ),
+      _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_GE_OQ));
+
+  // Compute approximate reciprocal sqrt.
+  Packet8f x = _mm256_rsqrt_ps(_x);
   // Do a single step of Newton's iteration.
-  x = pmul(x, pmadd(neg_half, pmul(x, x), p8f_one_point_five));
-
-  // Multiply the original _x by it's reciprocal square root to extract the
-  // square root.
-  return pmul(_x, x);
+  x = pmul(x, psub(pset1<Packet8f>(1.5f), pmul(half, pmul(x,x))));
+  // Flush results for denormals to zero.
+  return _mm256_andnot_ps(denormal_mask, pmul(_x,x));
 }
 #else
 template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h
index ac2fd8103..7b5f948e1 100644
--- a/Eigen/src/Core/arch/SSE/MathFunctions.h
+++ b/Eigen/src/Core/arch/SSE/MathFunctions.h
@@ -32,7 +32,7 @@ Packet4f plog<Packet4f>(const Packet4f& _x)
   /* the smallest non denormalized float number */
   _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos,  0x00800000);
   _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_inf,     0xff800000);//-1.f/0.f);
-  
+
   /* natural logarithm computed for 4 simultaneous float
     return NaN for x <= 0
   */
@@ -444,25 +444,33 @@ Packet4f pcos<Packet4f>(const Packet4f& _x)
 
 #if EIGEN_FAST_MATH
 
-// This is based on Quake3's fast inverse square root.
+// Functions for sqrt.
+// The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step
+// of Newton's method, at a cost of 1-2 bits of precision as opposed to the
+// exact solution. It does not handle +inf, or denormalized numbers correctly.
+// The main advantage of this approach is not just speed, but also the fact that
+// it can be inlined and pipelined with other computations, further reducing its
+// effective latency. This is similar to Quake3's fast inverse square root.
 // For detail see here: http://www.beyond3d.com/content/articles/8/
-// It lacks 1 (or 2 bits in some rare cases) of precision, and does not handle negative, +inf, or denormalized numbers correctly.
 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
 Packet4f psqrt<Packet4f>(const Packet4f& _x)
 {
   Packet4f half = pmul(_x, pset1<Packet4f>(.5f));
+  Packet4f denormal_mask = _mm_and_ps(
+      _mm_cmpge_ps(_x, _mm_setzero_ps()),
+      _mm_cmplt_ps(_x, pset1<Packet4f>((std::numeric_limits<float>::min)())));
 
-  /* select only the inverse sqrt of non-zero inputs */
-  Packet4f non_zero_mask = _mm_cmpge_ps(_x, pset1<Packet4f>((std::numeric_limits<float>::min)()));
-  Packet4f x = _mm_and_ps(non_zero_mask, _mm_rsqrt_ps(_x));
-
+  // Compute approximate reciprocal sqrt.
+  Packet4f x = _mm_rsqrt_ps(_x);
+  // Do a single step of Newton's iteration.
   x = pmul(x, psub(pset1<Packet4f>(1.5f), pmul(half, pmul(x,x))));
-  return pmul(_x,x);
+  // Flush results for denormals to zero.
+  return _mm_andnot_ps(denormal_mask, pmul(_x,x));
 }
 
 #else
 
-template<>EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 
+template<>EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
 Packet4f psqrt<Packet4f>(const Packet4f& x) { return _mm_sqrt_ps(x); }
 
 #endif
@@ -491,7 +499,7 @@ Packet4f prsqrt<Packet4f>(const Packet4f& _x) {
   Packet4f neg_mask = _mm_cmplt_ps(_x, _mm_setzero_ps());
   Packet4f zero_mask = _mm_andnot_ps(neg_mask, le_zero_mask);
   Packet4f infs_and_nans = _mm_or_ps(_mm_and_ps(neg_mask, p4f_nan),
-                                        _mm_and_ps(zero_mask, p4f_inf));
+                                     _mm_and_ps(zero_mask, p4f_inf));
 
   // Do a single step of Newton's iteration.
   x = pmul(x, pmadd(neg_half, pmul(x, x), p4f_one_point_five));
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index 1394d9f2b..c3d3e1521 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -193,7 +193,7 @@ template<typename Scalar> void packetmath()
     internal::pstore(data2+3*PacketSize, A3);
     VERIFY(areApprox(ref, data2, 4*PacketSize) && "internal::pbroadcast4");
   }
-  
+
   {
     for (int i=0; i<PacketSize*2; ++i)
       ref[i] = data1[i/PacketSize];
@@ -203,9 +203,9 @@ template<typename Scalar> void packetmath()
     internal::pstore(data2+1*PacketSize, A1);
     VERIFY(areApprox(ref, data2, 2*PacketSize) && "internal::pbroadcast2");
   }
- 
+
   VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst");
- 
+
   if(PacketSize>1)
   {
     for(int offset=0;offset<4;++offset)
@@ -315,7 +315,7 @@ template<typename Scalar> void packetmath_real()
   CHECK_CWISE1_IF(PacketTraits::HasRound, numext::round, internal::pround);
   CHECK_CWISE1_IF(PacketTraits::HasCeil, numext::ceil, internal::pceil);
   CHECK_CWISE1_IF(PacketTraits::HasFloor, numext::floor, internal::pfloor);
-  
+
   for (int i=0; i<size; ++i)
   {
     data1[i] = internal::random<Scalar>(-1,1);
@@ -440,12 +440,9 @@ template<typename Scalar> void packetmath_real()
     data1[0] = Scalar(-1.0f);
     h.store(data2, internal::plog(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
-#if !EIGEN_FAST_MATH
     h.store(data2, internal::psqrt(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
     VERIFY((numext::isnan)(data2[1]));
-#endif
-
   }
 }
 
@@ -459,7 +456,7 @@ template<typename Scalar> void packetmath_notcomplex()
   EIGEN_ALIGN_MAX Scalar data1[PacketTraits::size*4];
   EIGEN_ALIGN_MAX Scalar data2[PacketTraits::size*4];
   EIGEN_ALIGN_MAX Scalar ref[PacketTraits::size*4];
-  
+
   Array<Scalar,Dynamic,1>::Map(data1, PacketTraits::size*4).setRandom();
 
   ref[0] = data1[0];
@@ -478,7 +475,7 @@ template<typename Scalar> void packetmath_notcomplex()
   for (int i=0; i<PacketSize; ++i)
     ref[0] = (std::max)(ref[0],data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max");
-  
+
   for (int i=0; i<PacketSize; ++i)
     ref[i] = data1[0]+Scalar(i);
   internal::pstore(data2, internal::plset<Packet>(data1[0]));
@@ -490,12 +487,12 @@ template<typename Scalar,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar
   typedef internal::packet_traits<Scalar> PacketTraits;
   typedef typename PacketTraits::type Packet;
   const int PacketSize = PacketTraits::size;
-  
+
   internal::conj_if<ConjLhs> cj0;
   internal::conj_if<ConjRhs> cj1;
   internal::conj_helper<Scalar,Scalar,ConjLhs,ConjRhs> cj;
   internal::conj_helper<Packet,Packet,ConjLhs,ConjRhs> pcj;
-  
+
   for(int i=0;i<PacketSize;++i)
   {
     ref[i] = cj0(data1[i]) * cj1(data2[i]);
@@ -503,7 +500,7 @@ template<typename Scalar,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar
   }
   internal::pstore(pval,pcj.pmul(internal::pload<Packet>(data1),internal::pload<Packet>(data2)));
   VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmul");
-  
+
   for(int i=0;i<PacketSize;++i)
   {
     Scalar tmp = ref[i];
@@ -531,12 +528,12 @@ template<typename Scalar> void packetmath_complex()
     data1[i] = internal::random<Scalar>() * Scalar(1e2);
     data2[i] = internal::random<Scalar>() * Scalar(1e2);
   }
-  
+
   test_conj_helper<Scalar,false,false> (data1,data2,ref,pval);
   test_conj_helper<Scalar,false,true>  (data1,data2,ref,pval);
   test_conj_helper<Scalar,true,false>  (data1,data2,ref,pval);
   test_conj_helper<Scalar,true,true>   (data1,data2,ref,pval);
-  
+
   {
     for(int i=0;i<PacketSize;++i)
       ref[i] = Scalar(std::imag(data1[i]),std::real(data1[i]));
@@ -556,9 +553,9 @@ template<typename Scalar> void packetmath_scatter_gather()
   for (int i=0; i<PacketSize; ++i) {
     data1[i] = internal::random<Scalar>()/RealScalar(PacketSize);
   }
-  
+
   int stride = internal::random<int>(1,20);
-  
+
   EIGEN_ALIGN_MAX Scalar buffer[PacketSize*20];
   memset(buffer, 0, 20*sizeof(Packet));
   Packet packet = internal::pload<Packet>(data1);
@@ -594,7 +591,7 @@ void test_packetmath()
     CALL_SUBTEST_1( packetmath_notcomplex<float>() );
     CALL_SUBTEST_2( packetmath_notcomplex<double>() );
     CALL_SUBTEST_3( packetmath_notcomplex<int>() );
-    
+
     CALL_SUBTEST_1( packetmath_real<float>() );
     CALL_SUBTEST_2( packetmath_real<double>() );