Added special_packetmath test and tweaked bounds on tests.

Refactor shared packetmath code to header file.
(Squashed from PR !38)

4 files changed, 425 insertions, 313 deletions

diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index f81c07e40..578441f96 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -8,177 +8,7 @@
 // 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/.
 
-#include "main.h"
-#include "unsupported/Eigen/SpecialFunctions"
-#include <typeinfo>
-
-#if defined __GNUC__ && __GNUC__>=6
-  #pragma GCC diagnostic ignored "-Wignored-attributes"
-#endif
-// using namespace Eigen;
-
-#ifdef EIGEN_VECTORIZE_SSE
-const bool g_vectorize_sse = true;
-#else
-const bool g_vectorize_sse = false;
-#endif
-
-bool g_first_pass = true;
-
-namespace Eigen {
-namespace internal {
-
-template<typename T> T negate(const T& x) { return -x; }
-
-template<typename T>
-Map<const Array<unsigned char,sizeof(T),1> >
-bits(const T& x) {
-  return Map<const Array<unsigned char,sizeof(T),1> >(reinterpret_cast<const unsigned char *>(&x));
-}
-
-// The following implement bitwise operations on floating point types
-template<typename T,typename Bits,typename Func>
-T apply_bit_op(Bits a, Bits b, Func f) {
-  Array<unsigned char,sizeof(T),1> data;
-  T res;
-  for(Index i = 0; i < data.size(); ++i)
-    data[i] = f(a[i], b[i]);
-  // Note: The reinterpret_cast works around GCC's class-memaccess warnings:
-  std::memcpy(reinterpret_cast<unsigned char*>(&res), data.data(), sizeof(T));
-  return res;
-}
-
-#define EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,T)             \
-  template<> T EIGEN_CAT(p,OP)(const T& a,const T& b) { \
-    return apply_bit_op<T>(bits(a),bits(b),FUNC);     \
-  }
-
-#define EIGEN_TEST_MAKE_BITWISE(OP,FUNC)                  \
-  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,float)                 \
-  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,double)                \
-  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,half)                  \
-  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,std::complex<float>)   \
-  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,std::complex<double>)
-
-EIGEN_TEST_MAKE_BITWISE(xor,std::bit_xor<unsigned char>())
-EIGEN_TEST_MAKE_BITWISE(and,std::bit_and<unsigned char>())
-EIGEN_TEST_MAKE_BITWISE(or, std::bit_or<unsigned char>())
-struct bit_andnot{
-  template<typename T> T
-  operator()(T a, T b) const { return a & (~b); }
-};
-EIGEN_TEST_MAKE_BITWISE(andnot, bit_andnot())
-template<typename T>
-bool biteq(T a, T b) {
-  return (bits(a) == bits(b)).all();
-}
-
-}
-}
-
-// NOTE: we disable inlining for this function to workaround a GCC issue when using -O3 and the i387 FPU.
-template<typename Scalar> EIGEN_DONT_INLINE
-bool isApproxAbs(const Scalar& a, const Scalar& b, const typename NumTraits<Scalar>::Real& refvalue)
-{
-  return internal::isMuchSmallerThan(a-b, refvalue);
-}
-
-template<typename Scalar> bool areApproxAbs(const Scalar* a, const Scalar* b, int size, const typename NumTraits<Scalar>::Real& refvalue)
-{
-  for (int i=0; i<size; ++i)
-  {
-    if (!isApproxAbs(a[i],b[i],refvalue))
-    {
-      std::cout << "ref: [" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != vec: [" << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "]\n";
-      return false;
-    }
-  }
-  return true;
-}
-
-template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int size)
-{
-  for (int i=0; i<size; ++i)
-  {
-    if ((!internal::biteq(a[i],b[i])) && a[i]!=b[i] && !internal::isApprox(a[i],b[i]))
-    {
-      std::cout << "ref: [" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != vec: [" << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "]\n";
-      return false;
-    }
-  }
-  return true;
-}
-
-#define CHECK_CWISE1(REFOP, POP) { \
-  for (int i=0; i<PacketSize; ++i) \
-    ref[i] = REFOP(data1[i]); \
-  internal::pstore(data2, POP(internal::pload<Packet>(data1))); \
-  VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
-}
-
-template<bool Cond,typename Packet>
-struct packet_helper
-{
-  template<typename T>
-  inline Packet load(const T* from) const { return internal::pload<Packet>(from); }
-
-  template<typename T>
-  inline Packet loadu(const T* from) const { return internal::ploadu<Packet>(from); }
-
-  template<typename T>
-  inline Packet load(const T* from, unsigned long long umask) const { return internal::ploadu<Packet>(from, umask); }
-
-  template<typename T>
-  inline void store(T* to, const Packet& x) const { internal::pstore(to,x); }
-
-  template<typename T>
-  inline void store(T* to, const Packet& x, unsigned long long umask) const { internal::pstoreu(to, x, umask); }
-};
-
-template<typename Packet>
-struct packet_helper<false,Packet>
-{
-  template<typename T>
-  inline T load(const T* from) const { return *from; }
-
-  template<typename T>
-  inline T loadu(const T* from) const { return *from; }
-
-  template<typename T>
-  inline T load(const T* from, unsigned long long) const { return *from; }
-
-  template<typename T>
-  inline void store(T* to, const T& x) const { *to = x; }
-
-  template<typename T>
-  inline void store(T* to, const T& x, unsigned long long) const { *to = x; }
-};
-
-#define CHECK_CWISE1_IF(COND, REFOP, POP) if(COND) { \
-  packet_helper<COND,Packet> h; \
-  for (int i=0; i<PacketSize; ++i) \
-    ref[i] = REFOP(data1[i]); \
-  h.store(data2, POP(h.load(data1))); \
-  VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
-}
-
-#define CHECK_CWISE2_IF(COND, REFOP, POP) if(COND) { \
-  packet_helper<COND,Packet> h; \
-  for (int i=0; i<PacketSize; ++i) \
-    ref[i] = REFOP(data1[i], data1[i+PacketSize]); \
-  h.store(data2, POP(h.load(data1),h.load(data1+PacketSize))); \
-  VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
-}
-
-#define CHECK_CWISE3_IF(COND, REFOP, POP) if (COND) {                      \
-  packet_helper<COND, Packet> h;                                           \
-  for (int i = 0; i < PacketSize; ++i)                                     \
-    ref[i] =                                                               \
-        REFOP(data1[i], data1[i + PacketSize], data1[i + 2 * PacketSize]); \
-  h.store(data2, POP(h.load(data1), h.load(data1 + PacketSize),            \
-                     h.load(data1 + 2 * PacketSize)));                     \
-  VERIFY(areApprox(ref, data2, PacketSize) && #POP);                       \
-}
+#include "packetmath_test_shared.h"
 
 #define REF_ADD(a,b) ((a)+(b))
 #define REF_SUB(a,b) ((a)-(b))
@@ -213,23 +43,23 @@ template<typename Scalar,typename Packet> void packetmath()
   }
 
   internal::pstore(data2, internal::pload<Packet>(data1));
-  VERIFY(areApprox(data1, data2, PacketSize) && "aligned load/store");
+  VERIFY(test::areApprox(data1, data2, PacketSize) && "aligned load/store");
 
   for (int offset=0; offset<PacketSize; ++offset)
   {
     internal::pstore(data2, internal::ploadu<Packet>(data1+offset));
-    VERIFY(areApprox(data1+offset, data2, PacketSize) && "internal::ploadu");
+    VERIFY(test::areApprox(data1+offset, data2, PacketSize) && "internal::ploadu");
   }
 
   for (int offset=0; offset<PacketSize; ++offset)
   {
     internal::pstoreu(data2+offset, internal::pload<Packet>(data1));
-    VERIFY(areApprox(data1, data2+offset, PacketSize) && "internal::pstoreu");
+    VERIFY(test::areApprox(data1, data2+offset, PacketSize) && "internal::pstoreu");
   }
 
   if (internal::unpacket_traits<Packet>::masked_load_available)
   {
-    packet_helper<internal::unpacket_traits<Packet>::masked_load_available, Packet> h;
+    test::packet_helper<internal::unpacket_traits<Packet>::masked_load_available, Packet> h;
     unsigned long long max_umask = (0x1ull << PacketSize);
 
     for (int offset=0; offset<PacketSize; ++offset)
@@ -239,14 +69,14 @@ template<typename Scalar,typename Packet> void packetmath()
         h.store(data2, h.load(data1+offset, umask));
         for (int k=0; k<PacketSize; ++k)
           data3[k] = ((umask & ( 0x1ull << k )) >> k) ? data1[k+offset] : Scalar(0);
-        VERIFY(areApprox(data3, data2, PacketSize) && "internal::ploadu masked");
+        VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::ploadu masked");
       }
     }
   }
 
   if (internal::unpacket_traits<Packet>::masked_store_available)
   {
-    packet_helper<internal::unpacket_traits<Packet>::masked_store_available, Packet> h;
+    test::packet_helper<internal::unpacket_traits<Packet>::masked_store_available, Packet> h;
     unsigned long long max_umask = (0x1ull << PacketSize);
 
     for (int offset=0; offset<PacketSize; ++offset)
@@ -257,7 +87,7 @@ template<typename Scalar,typename Packet> void packetmath()
         h.store(data2, h.loadu(data1+offset), umask);
         for (int k=0; k<PacketSize; ++k)
           data3[k] = ((umask & ( 0x1ull << k )) >> k) ? data1[k+offset] : Scalar(0);
-        VERIFY(areApprox(data3, data2, PacketSize) && "internal::pstoreu masked");
+        VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::pstoreu masked");
       }
     }
   }
@@ -290,7 +120,7 @@ template<typename Scalar,typename Packet> void packetmath()
 
     // palign is not used anymore, so let's just put a warning if it fails
     ++g_test_level;
-    VERIFY(areApprox(ref, data2, PacketSize) && "internal::palign");
+    VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::palign");
     --g_test_level;
   }
 
@@ -317,7 +147,7 @@ template<typename Scalar,typename Packet> void packetmath()
     for (int i=0; i<PacketSize; ++i)
       ref[i] = data1[offset];
     internal::pstore(data2, internal::pset1<Packet>(data1[offset]));
-    VERIFY(areApprox(ref, data2, PacketSize) && "internal::pset1");
+    VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::pset1");
   }
 
   {
@@ -329,7 +159,7 @@ template<typename Scalar,typename Packet> void packetmath()
     internal::pstore(data2+1*PacketSize, A1);
     internal::pstore(data2+2*PacketSize, A2);
     internal::pstore(data2+3*PacketSize, A3);
-    VERIFY(areApprox(ref, data2, 4*PacketSize) && "internal::pbroadcast4");
+    VERIFY(test::areApprox(ref, data2, 4*PacketSize) && "internal::pbroadcast4");
   }
 
   {
@@ -339,7 +169,7 @@ template<typename Scalar,typename Packet> void packetmath()
     internal::pbroadcast2<Packet>(data1, A0, A1);
     internal::pstore(data2+0*PacketSize, A0);
     internal::pstore(data2+1*PacketSize, A1);
-    VERIFY(areApprox(ref, data2, 2*PacketSize) && "internal::pbroadcast2");
+    VERIFY(test::areApprox(ref, data2, 2*PacketSize) && "internal::pbroadcast2");
   }
 
   VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst");
@@ -352,7 +182,7 @@ template<typename Scalar,typename Packet> void packetmath()
       for(int i=0;i<PacketSize/2;++i)
         ref[2*i+0] = ref[2*i+1] = data1[offset+i];
       internal::pstore(data2,internal::ploaddup<Packet>(data1+offset));
-      VERIFY(areApprox(ref, data2, PacketSize) && "ploaddup");
+      VERIFY(test::areApprox(ref, data2, PacketSize) && "ploaddup");
     }
   }
 
@@ -364,14 +194,14 @@ template<typename Scalar,typename Packet> void packetmath()
       for(int i=0;i<PacketSize/4;++i)
         ref[4*i+0] = ref[4*i+1] = ref[4*i+2] = ref[4*i+3] = data1[offset+i];
       internal::pstore(data2,internal::ploadquad<Packet>(data1+offset));
-      VERIFY(areApprox(ref, data2, PacketSize) && "ploadquad");
+      VERIFY(test::areApprox(ref, data2, PacketSize) && "ploadquad");
     }
   }
 
   ref[0] = Scalar(0);
   for (int i=0; i<PacketSize; ++i)
     ref[0] += data1[i];
-  VERIFY(isApproxAbs(ref[0], internal::predux(internal::pload<Packet>(data1)), refvalue) && "internal::predux");
+  VERIFY(test::isApproxAbs(ref[0], internal::predux(internal::pload<Packet>(data1)), refvalue) && "internal::predux");
 
   if(PacketSize==8 && internal::unpacket_traits<typename internal::unpacket_traits<Packet>::half>::size ==4) // so far, predux_half_downto4 is only required in such a case
   {
@@ -381,7 +211,7 @@ template<typename Scalar,typename Packet> void packetmath()
     for (int i=0; i<PacketSize; ++i)
       ref[i%HalfPacketSize] += data1[i];
     internal::pstore(data2, internal::predux_half_dowto4(internal::pload<Packet>(data1)));
-    VERIFY(areApprox(ref, data2, HalfPacketSize) && "internal::predux_half_dowto4");
+    VERIFY(test::areApprox(ref, data2, HalfPacketSize) && "internal::predux_half_dowto4");
   }
 
   ref[0] = Scalar(1);
@@ -399,13 +229,13 @@ template<typename Scalar,typename Packet> void packetmath()
       packets[j] = internal::pload<Packet>(data1+j*PacketSize);
     }
     internal::pstore(data2, internal::preduxp(packets));
-    VERIFY(areApproxAbs(ref, data2, PacketSize, refvalue) && "internal::preduxp");
+    VERIFY(test::areApproxAbs(ref, data2, PacketSize, refvalue) && "internal::preduxp");
   }
 
   for (int i=0; i<PacketSize; ++i)
     ref[i] = data1[PacketSize-i-1];
   internal::pstore(data2, internal::preverse(internal::pload<Packet>(data1)));
-  VERIFY(areApprox(ref, data2, PacketSize) && "internal::preverse");
+  VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::preverse");
 
   internal::PacketBlock<Packet> kernel;
   for (int i=0; i<PacketSize; ++i) {
@@ -415,7 +245,7 @@ template<typename Scalar,typename Packet> void packetmath()
   for (int i=0; i<PacketSize; ++i) {
     internal::pstore(data2, kernel.packet[i]);
     for (int j = 0; j < PacketSize; ++j) {
-      VERIFY(isApproxAbs(data2[j], data1[i+j*PacketSize], refvalue) && "ptranspose");
+      VERIFY(test::isApproxAbs(data2[j], data1[i+j*PacketSize], refvalue) && "ptranspose");
     }
   }
 
@@ -431,7 +261,7 @@ template<typename Scalar,typename Packet> void packetmath()
     EIGEN_ALIGN_MAX Scalar result[size];
     internal::pstore(result, blend);
     for (int i = 0; i < PacketSize; ++i) {
-      VERIFY(isApproxAbs(result[i], (selector.select[i] ? data1[i] : data2[i]), refvalue));
+      VERIFY(test::isApproxAbs(result[i], (selector.select[i] ? data1[i] : data2[i]), refvalue));
     }
   }
 
@@ -442,7 +272,7 @@ template<typename Scalar,typename Packet> void packetmath()
     Scalar s = internal::random<Scalar>();
     ref[0] = s;
     internal::pstore(data2, internal::pinsertfirst(internal::pload<Packet>(data1),s));
-    VERIFY(areApprox(ref, data2, PacketSize) && "internal::pinsertfirst");
+    VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::pinsertfirst");
   }
 
   if (PacketTraits::HasBlend || g_vectorize_sse) {
@@ -452,7 +282,7 @@ template<typename Scalar,typename Packet> void packetmath()
     Scalar s = internal::random<Scalar>();
     ref[PacketSize-1] = s;
     internal::pstore(data2, internal::pinsertlast(internal::pload<Packet>(data1),s));
-    VERIFY(areApprox(ref, data2, PacketSize) && "internal::pinsertlast");
+    VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::pinsertlast");
   }
 
   {
@@ -508,6 +338,19 @@ template<typename Scalar,typename Packet> void packetmath_real()
 
   for (int i=0; i<size; ++i)
   {
+    data1[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
+    data2[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
+  }
+
+  if(internal::random<float>(0,1)<0.1f)
+     data1[internal::random<int>(0, PacketSize)] = 0;
+
+  CHECK_CWISE1_IF(PacketTraits::HasSqrt, std::sqrt, internal::psqrt);
+  CHECK_CWISE1_IF(PacketTraits::HasLog, std::log, internal::plog);
+  CHECK_CWISE1_IF(PacketTraits::HasRsqrt, Scalar(1)/std::sqrt, internal::prsqrt);
+
+  for (int i=0; i<size; ++i)
+  {
     data1[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-3,3));
     data2[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-3,3));
   }
@@ -554,7 +397,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
   {
     data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
     data1[1] = std::numeric_limits<Scalar>::epsilon();
-    packet_helper<PacketTraits::HasExp,Packet> h;
+    test::packet_helper<PacketTraits::HasExp,Packet> h;
     h.store(data2, internal::pexp(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
     VERIFY_IS_EQUAL(std::exp(std::numeric_limits<Scalar>::epsilon()), data2[1]);
@@ -581,77 +424,12 @@ template<typename Scalar,typename Packet> void packetmath_real()
   if (PacketTraits::HasTanh) {
     // NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details.
     data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
-    packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h;
+    test::packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h;
     h.store(data2, internal::ptanh(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
   }
 
-#if EIGEN_HAS_C99_MATH
-  {
-    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
-    packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;
-    h.store(data2, internal::plgamma(h.load(data1)));
-    VERIFY((numext::isnan)(data2[0]));
-  }
-  if (internal::packet_traits<Scalar>::HasErf) {
-    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
-    packet_helper<internal::packet_traits<Scalar>::HasErf,Packet> h;
-    h.store(data2, internal::perf(h.load(data1)));
-    VERIFY((numext::isnan)(data2[0]));
-  }
-  {
-    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
-    packet_helper<internal::packet_traits<Scalar>::HasErfc,Packet> h;
-    h.store(data2, internal::perfc(h.load(data1)));
-    VERIFY((numext::isnan)(data2[0]));
-  }
-  {
-    for (int i=0; i<size; ++i) {
-      data1[i] = internal::random<Scalar>(0,1);
-    }
-    CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasNdtri, numext::ndtri, internal::pndtri);
-  }
-#endif  // EIGEN_HAS_C99_MATH
-
-  for (int i=0; i<size; ++i)
-  {
-    data1[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
-    data2[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
-  }
-
-  if(internal::random<float>(0,1)<0.1f)
-     data1[internal::random<int>(0, PacketSize)] = 0;
-  CHECK_CWISE1_IF(PacketTraits::HasSqrt, std::sqrt, internal::psqrt);
-  CHECK_CWISE1_IF(PacketTraits::HasLog, std::log, internal::plog);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1);
-
-  data1[0] = std::numeric_limits<Scalar>::infinity();
-  CHECK_CWISE1_IF(PacketTraits::HasRsqrt, Scalar(1)/std::sqrt, internal::prsqrt);
-
-  // Use a smaller data range for the positive bessel operations as these
-  // can have much more error at very small and very large values.
-  for (int i=0; i<size; ++i) {
-      data1[i] = internal::random<Scalar>(0.01,1) * std::pow(
-          Scalar(10), internal::random<Scalar>(-1,2));
-      data2[i] = internal::random<Scalar>(0.01,1) * std::pow(
-          Scalar(10), internal::random<Scalar>(-1,2));
-  }
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y0, internal::pbessel_y0);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y1, internal::pbessel_y1);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0, internal::pbessel_k0);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0e, internal::pbessel_k0e);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1, internal::pbessel_k1);
-  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1e, internal::pbessel_k1e);
-
 #if EIGEN_HAS_C99_MATH && (__cplusplus > 199711L)
-  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
-  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
-  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc);
   data1[0] = std::numeric_limits<Scalar>::infinity();
   data1[1] = Scalar(-1);
   CHECK_CWISE1_IF(PacketTraits::HasLog1p, std::log1p, internal::plog1p);
@@ -666,7 +444,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
     data1[1] = std::numeric_limits<Scalar>::epsilon();
     if(PacketTraits::HasLog)
     {
-      packet_helper<PacketTraits::HasLog,Packet> h;
+      test::packet_helper<PacketTraits::HasLog,Packet> h;
       h.store(data2, internal::plog(h.load(data1)));
       VERIFY((numext::isnan)(data2[0]));
       VERIFY_IS_EQUAL(std::log(std::numeric_limits<Scalar>::epsilon()), data2[1]);
@@ -698,7 +476,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
       VERIFY((numext::isinf)(data2[0]));
     }
     if(PacketTraits::HasLog1p) {
-      packet_helper<PacketTraits::HasLog1p,Packet> h;
+      test::packet_helper<PacketTraits::HasLog1p,Packet> h;
       data1[0] = Scalar(-2);
       data1[1] = -std::numeric_limits<Scalar>::infinity();
       h.store(data2, internal::plog1p(h.load(data1)));
@@ -707,7 +485,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
     }
     if(PacketTraits::HasSqrt)
     {
-      packet_helper<PacketTraits::HasSqrt,Packet> h;
+      test::packet_helper<PacketTraits::HasSqrt,Packet> h;
       data1[0] = Scalar(-1.0f);
       data1[1] = -std::numeric_limits<Scalar>::denorm_min();
       h.store(data2, internal::psqrt(h.load(data1)));
@@ -716,7 +494,7 @@ template<typename Scalar,typename Packet> void packetmath_real()
     }
     if(PacketTraits::HasCos)
     {
-      packet_helper<PacketTraits::HasCos,Packet> h;
+      test::packet_helper<PacketTraits::HasCos,Packet> h;
       for(Scalar k = 1; k<Scalar(10000)/std::numeric_limits<Scalar>::epsilon(); k*=2)
       {
         for(int k1=0;k1<=1; ++k1)
@@ -792,7 +570,7 @@ template<typename Scalar,typename Packet> void packetmath_notcomplex()
   for (int i=0; i<PacketSize; ++i)
     ref[i] = data1[0]+Scalar(i);
   internal::pstore(data2, internal::plset<Packet>(data1[0]));
-  VERIFY(areApprox(ref, data2, PacketSize) && "internal::plset");
+  VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::plset");
 
   {
     unsigned char* data1_bits = reinterpret_cast<unsigned char*>(data1);
@@ -833,7 +611,7 @@ template<typename Scalar,typename Packet,bool ConjLhs,bool ConjRhs> void test_co
     VERIFY(internal::isApprox(ref[i], cj.pmul(data1[i],data2[i])) && "conj_helper pmul");
   }
   internal::pstore(pval,pcj.pmul(internal::pload<Packet>(data1),internal::pload<Packet>(data2)));
-  VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmul");
+  VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmul");
 
   for(int i=0;i<PacketSize;++i)
   {
@@ -842,7 +620,7 @@ template<typename Scalar,typename Packet,bool ConjLhs,bool ConjRhs> void test_co
     VERIFY(internal::isApprox(ref[i], cj.pmadd(data1[i],data2[i],tmp)) && "conj_helper pmadd");
   }
   internal::pstore(pval,pcj.pmadd(internal::pload<Packet>(data1),internal::pload<Packet>(data2),internal::pload<Packet>(pval)));
-  VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmadd");
+  VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmadd");
 }
 
 template<typename Scalar,typename Packet> void packetmath_complex()
@@ -870,7 +648,7 @@ template<typename Scalar,typename Packet> void packetmath_complex()
     for(int i=0;i<PacketSize;++i)
       ref[i] = Scalar(std::imag(data1[i]),std::real(data1[i]));
     internal::pstore(pval,internal::pcplxflip(internal::pload<Packet>(data1)));
-    VERIFY(areApprox(ref, pval, PacketSize) && "pcplxflip");
+    VERIFY(test::areApprox(ref, pval, PacketSize) && "pcplxflip");
   }
 }
 
@@ -893,9 +671,11 @@ template<typename Scalar,typename Packet> void packetmath_scatter_gather()
 
   for (int i = 0; i < PacketSize*20; ++i) {
     if ((i%stride) == 0 && i<stride*PacketSize) {
-      VERIFY(isApproxAbs(buffer[i], data1[i/stride], refvalue) && "pscatter");
+      VERIFY(
+          test::isApproxAbs(buffer[i], data1[i/stride], refvalue) && "pscatter");
     } else {
-      VERIFY(isApproxAbs(buffer[i], Scalar(0), refvalue) && "pscatter");
+      VERIFY(
+          test::isApproxAbs(buffer[i], Scalar(0), refvalue) && "pscatter");
     }
   }
 
@@ -905,17 +685,12 @@ template<typename Scalar,typename Packet> void packetmath_scatter_gather()
   packet = internal::pgather<Scalar, Packet>(buffer, 7);
   internal::pstore(data1, packet);
   for (int i = 0; i < PacketSize; ++i) {
-    VERIFY(isApproxAbs(data1[i], buffer[i*7], refvalue) && "pgather");
+    VERIFY(test::isApproxAbs(data1[i], buffer[i*7], refvalue) && "pgather");
   }
 }
 
-
-template<
-  typename Scalar,
-  typename PacketType,
-  bool IsComplex = NumTraits<Scalar>::IsComplex,
-  bool IsInteger = NumTraits<Scalar>::IsInteger>
-struct runall;
+namespace Eigen {
+namespace test {
 
 template<typename Scalar,typename PacketType>
 struct runall<Scalar,PacketType,false,false> { // i.e. float or double
@@ -945,49 +720,20 @@ struct runall<Scalar,PacketType,true,false> { // i.e. complex
   }
 };
 
-template<
-  typename Scalar,
-  typename PacketType = typename internal::packet_traits<Scalar>::type,
-  bool Vectorized = internal::packet_traits<Scalar>::Vectorizable,
-  bool HasHalf = !internal::is_same<typename internal::unpacket_traits<PacketType>::half,PacketType>::value >
-struct runner;
-
-template<typename Scalar,typename PacketType>
-struct runner<Scalar,PacketType,true,true>
-{
-  static void run() {
-    runall<Scalar,PacketType>::run();
-    runner<Scalar,typename internal::unpacket_traits<PacketType>::half>::run();
-  }
-};
-
-template<typename Scalar,typename PacketType>
-struct runner<Scalar,PacketType,true,false>
-{
-  static void run() {
-    runall<Scalar,PacketType>::run();
-    runall<Scalar,Scalar>::run();
-  }
-};
+}
+}
 
-template<typename Scalar,typename PacketType>
-struct runner<Scalar,PacketType,false,false>
-{
-  static void run() {
-    runall<Scalar,PacketType>::run();
-  }
-};
 
 EIGEN_DECLARE_TEST(packetmath)
 {
   g_first_pass = true;
   for(int i = 0; i < g_repeat; i++) {
 
-    CALL_SUBTEST_1( runner<float>::run() );
-    CALL_SUBTEST_2( runner<double>::run() );
-    CALL_SUBTEST_3( runner<int>::run() );
-    CALL_SUBTEST_4( runner<std::complex<float> >::run() );
-    CALL_SUBTEST_5( runner<std::complex<double> >::run() );
+    CALL_SUBTEST_1( test::runner<float>::run() );
+    CALL_SUBTEST_2( test::runner<double>::run() );
+    CALL_SUBTEST_3( test::runner<int>::run() );
+    CALL_SUBTEST_4( test::runner<std::complex<float> >::run() );
+    CALL_SUBTEST_5( test::runner<std::complex<double> >::run() );
     CALL_SUBTEST_6(( packetmath<half,internal::packet_traits<half>::type>() ));
     g_first_pass = false;
   }
diff --git a/test/packetmath_test_shared.h b/test/packetmath_test_shared.h
new file mode 100644
index 000000000..046fd8104
--- /dev/null
+++ b/test/packetmath_test_shared.h
@@ -0,0 +1,225 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// 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/.
+
+#include "main.h"
+#include <typeinfo>
+
+#if defined __GNUC__ && __GNUC__>=6
+  #pragma GCC diagnostic ignored "-Wignored-attributes"
+#endif
+// using namespace Eigen;
+
+#ifdef EIGEN_VECTORIZE_SSE
+const bool g_vectorize_sse = true;
+#else
+const bool g_vectorize_sse = false;
+#endif
+
+bool g_first_pass = true;
+
+namespace Eigen {
+namespace internal {
+
+template<typename T> T negate(const T& x) { return -x; }
+
+template<typename T>
+Map<const Array<unsigned char,sizeof(T),1> >
+bits(const T& x) {
+  return Map<const Array<unsigned char,sizeof(T),1> >(reinterpret_cast<const unsigned char *>(&x));
+}
+
+// The following implement bitwise operations on floating point types
+template<typename T,typename Bits,typename Func>
+T apply_bit_op(Bits a, Bits b, Func f) {
+  Array<unsigned char,sizeof(T),1> data;
+  T res;
+  for(Index i = 0; i < data.size(); ++i)
+    data[i] = f(a[i], b[i]);
+  // Note: The reinterpret_cast works around GCC's class-memaccess warnings:
+  std::memcpy(reinterpret_cast<unsigned char*>(&res), data.data(), sizeof(T));
+  return res;
+}
+
+#define EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,T)             \
+  template<> T EIGEN_CAT(p,OP)(const T& a,const T& b) { \
+    return apply_bit_op<T>(bits(a),bits(b),FUNC);     \
+  }
+
+#define EIGEN_TEST_MAKE_BITWISE(OP,FUNC)                  \
+  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,float)                 \
+  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,double)                \
+  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,half)                  \
+  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,std::complex<float>)   \
+  EIGEN_TEST_MAKE_BITWISE2(OP,FUNC,std::complex<double>)
+
+EIGEN_TEST_MAKE_BITWISE(xor,std::bit_xor<unsigned char>())
+EIGEN_TEST_MAKE_BITWISE(and,std::bit_and<unsigned char>())
+EIGEN_TEST_MAKE_BITWISE(or, std::bit_or<unsigned char>())
+struct bit_andnot{
+  template<typename T> T
+  operator()(T a, T b) const { return a & (~b); }
+};
+EIGEN_TEST_MAKE_BITWISE(andnot, bit_andnot())
+template<typename T>
+bool biteq(T a, T b) {
+  return (bits(a) == bits(b)).all();
+}
+
+}
+
+namespace test {
+
+// NOTE: we disable inlining for this function to workaround a GCC issue when using -O3 and the i387 FPU.
+template<typename Scalar> EIGEN_DONT_INLINE
+bool isApproxAbs(const Scalar& a, const Scalar& b, const typename NumTraits<Scalar>::Real& refvalue)
+{
+  return internal::isMuchSmallerThan(a-b, refvalue);
+}
+
+template<typename Scalar> bool areApproxAbs(const Scalar* a, const Scalar* b, int size, const typename NumTraits<Scalar>::Real& refvalue)
+{
+  for (int i=0; i<size; ++i)
+  {
+    if (!isApproxAbs(a[i],b[i],refvalue))
+    {
+      std::cout << "ref: [" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != vec: [" << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "]\n";
+      return false;
+    }
+  }
+  return true;
+}
+
+template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int size)
+{
+  for (int i=0; i<size; ++i)
+  {
+    if ((!internal::biteq(a[i],b[i])) && a[i]!=b[i] && !internal::isApprox(a[i],b[i]))
+    {
+      std::cout << "ref: [" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != vec: [" << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "]\n";
+      return false;
+    }
+  }
+  return true;
+}
+
+#define CHECK_CWISE1(REFOP, POP) { \
+  for (int i=0; i<PacketSize; ++i) \
+    ref[i] = REFOP(data1[i]); \
+  internal::pstore(data2, POP(internal::pload<Packet>(data1))); \
+  VERIFY(test::areApprox(ref, data2, PacketSize) && #POP); \
+}
+
+template<bool Cond,typename Packet>
+struct packet_helper
+{
+  template<typename T>
+  inline Packet load(const T* from) const { return internal::pload<Packet>(from); }
+
+  template<typename T>
+  inline Packet loadu(const T* from) const { return internal::ploadu<Packet>(from); }
+
+  template<typename T>
+  inline Packet load(const T* from, unsigned long long umask) const { return internal::ploadu<Packet>(from, umask); }
+
+  template<typename T>
+  inline void store(T* to, const Packet& x) const { internal::pstore(to,x); }
+
+  template<typename T>
+  inline void store(T* to, const Packet& x, unsigned long long umask) const { internal::pstoreu(to, x, umask); }
+};
+
+template<typename Packet>
+struct packet_helper<false,Packet>
+{
+  template<typename T>
+  inline T load(const T* from) const { return *from; }
+
+  template<typename T>
+  inline T loadu(const T* from) const { return *from; }
+
+  template<typename T>
+  inline T load(const T* from, unsigned long long) const { return *from; }
+
+  template<typename T>
+  inline void store(T* to, const T& x) const { *to = x; }
+
+  template<typename T>
+  inline void store(T* to, const T& x, unsigned long long) const { *to = x; }
+};
+
+#define CHECK_CWISE1_IF(COND, REFOP, POP) if(COND) { \
+  test::packet_helper<COND,Packet> h; \
+  for (int i=0; i<PacketSize; ++i) \
+    ref[i] = REFOP(data1[i]); \
+  h.store(data2, POP(h.load(data1))); \
+  VERIFY(test::areApprox(ref, data2, PacketSize) && #POP); \
+}
+
+#define CHECK_CWISE2_IF(COND, REFOP, POP) if(COND) { \
+  test::packet_helper<COND,Packet> h; \
+  for (int i=0; i<PacketSize; ++i) \
+    ref[i] = REFOP(data1[i], data1[i+PacketSize]); \
+  h.store(data2, POP(h.load(data1),h.load(data1+PacketSize))); \
+  VERIFY(test::areApprox(ref, data2, PacketSize) && #POP); \
+}
+
+#define CHECK_CWISE3_IF(COND, REFOP, POP) if (COND) {                      \
+  test::packet_helper<COND, Packet> h;                                     \
+  for (int i = 0; i < PacketSize; ++i)                                     \
+    ref[i] =                                                               \
+        REFOP(data1[i], data1[i + PacketSize], data1[i + 2 * PacketSize]); \
+  h.store(data2, POP(h.load(data1), h.load(data1 + PacketSize),            \
+                     h.load(data1 + 2 * PacketSize)));                     \
+  VERIFY(test::areApprox(ref, data2, PacketSize) && #POP);                 \
+}
+
+// Specialize the runall struct in your test file by defining run().
+template<
+  typename Scalar,
+  typename PacketType,
+  bool IsComplex = NumTraits<Scalar>::IsComplex,
+  bool IsInteger = NumTraits<Scalar>::IsInteger>
+struct runall;
+
+template<
+  typename Scalar,
+  typename PacketType = typename internal::packet_traits<Scalar>::type,
+  bool Vectorized = internal::packet_traits<Scalar>::Vectorizable,
+  bool HasHalf = !internal::is_same<typename internal::unpacket_traits<PacketType>::half,PacketType>::value >
+struct runner;
+
+template<typename Scalar,typename PacketType>
+struct runner<Scalar,PacketType,true,true>
+{
+  static void run() {
+    runall<Scalar,PacketType>::run();
+    runner<Scalar,typename internal::unpacket_traits<PacketType>::half>::run();
+  }
+};
+
+template<typename Scalar,typename PacketType>
+struct runner<Scalar,PacketType,true,false>
+{
+  static void run() {
+    runall<Scalar,PacketType>::run();
+    runall<Scalar,Scalar>::run();
+  }
+};
+
+template<typename Scalar,typename PacketType>
+struct runner<Scalar,PacketType,false,false>
+{
+  static void run() {
+    runall<Scalar,PacketType>::run();
+  }
+};
+
+}
+}
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 9db965ad8..298db04ea 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -108,6 +108,7 @@ ei_add_test(levenberg_marquardt)
 ei_add_test(kronecker_product)
 ei_add_test(bessel_functions)
 ei_add_test(special_functions)
+ei_add_test(special_packetmath "-DEIGEN_FAST_MATH=1")
 
 if(EIGEN_TEST_CXX11)
   if(EIGEN_TEST_SYCL)
diff --git a/unsupported/test/special_packetmath.cpp b/unsupported/test/special_packetmath.cpp
new file mode 100644
index 000000000..87b87351b
--- /dev/null
+++ b/unsupported/test/special_packetmath.cpp
@@ -0,0 +1,140 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// 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/.
+
+#include "packetmath_test_shared.h"
+#include "../Eigen/SpecialFunctions"
+
+template<typename Scalar,typename Packet> void packetmath_real()
+{
+  using std::abs;
+  typedef internal::packet_traits<Scalar> PacketTraits;
+  const int PacketSize = internal::unpacket_traits<Packet>::size;
+
+  const int size = PacketSize*4;
+  EIGEN_ALIGN_MAX Scalar data1[PacketSize*4];
+  EIGEN_ALIGN_MAX Scalar data2[PacketSize*4];
+  EIGEN_ALIGN_MAX Scalar ref[PacketSize*4];
+
+#if EIGEN_HAS_C99_MATH
+  {
+    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
+    test::packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;
+    h.store(data2, internal::plgamma(h.load(data1)));
+    VERIFY((numext::isnan)(data2[0]));
+  }
+  if (internal::packet_traits<Scalar>::HasErf) {
+    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
+    test::packet_helper<internal::packet_traits<Scalar>::HasErf,Packet> h;
+    h.store(data2, internal::perf(h.load(data1)));
+    VERIFY((numext::isnan)(data2[0]));
+  }
+  {
+    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
+    test::packet_helper<internal::packet_traits<Scalar>::HasErfc,Packet> h;
+    h.store(data2, internal::perfc(h.load(data1)));
+    VERIFY((numext::isnan)(data2[0]));
+  }
+  {
+    for (int i=0; i<size; ++i) {
+      data1[i] = internal::random<Scalar>(0,1);
+    }
+    CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasNdtri, numext::ndtri, internal::pndtri);
+  }
+#endif  // EIGEN_HAS_C99_MATH
+
+  // For bessel_i*e and bessel_j*, the valid range is negative reals.
+  for (int i=0; i<size; ++i)
+  {
+    data1[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
+    data2[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6));
+  }
+
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1);
+
+  // Use a smaller data range for the bessel_i* as these can become very large.
+  // Following #1693, we also restrict this range further to avoid inf's due to
+  // differences in pexp and exp.
+  for (int i=0; i<size; ++i) {
+      data1[i] = internal::random<Scalar>(0.01,1) * std::pow(
+          Scalar(9), internal::random<Scalar>(-1,2));
+      data2[i] = internal::random<Scalar>(0.01,1) * std::pow(
+          Scalar(9), internal::random<Scalar>(-1,2));
+  }
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1);
+
+
+  // y_i, and k_i are valid for x > 0.
+  for (int i=0; i<size; ++i)
+  {
+    data1[i] = internal::random<Scalar>(0.01,1) * std::pow(Scalar(10), internal::random<Scalar>(-2,5));
+    data2[i] = internal::random<Scalar>(0.01,1) * std::pow(Scalar(10), internal::random<Scalar>(-2,5));
+  }
+
+  // TODO(srvasude): Re-enable this test once properly investigated why the
+  // scalar and vector paths differ.
+  // CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y0, internal::pbessel_y0);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y1, internal::pbessel_y1);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0e, internal::pbessel_k0e);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1e, internal::pbessel_k1e);
+
+  // Following #1693, we restrict the range for exp to avoid zeroing out too
+  // fast.
+  for (int i=0; i<size; ++i) {
+      data1[i] = internal::random<Scalar>(0.01,1) * std::pow(
+          Scalar(9), internal::random<Scalar>(-1,2));
+      data2[i] = internal::random<Scalar>(0.01,1) * std::pow(
+          Scalar(9), internal::random<Scalar>(-1,2));
+  }
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0, internal::pbessel_k0);
+  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1, internal::pbessel_k1);
+
+
+  for (int i=0; i<size; ++i) {
+      data1[i] = internal::random<Scalar>(0.01,1) * std::pow(
+          Scalar(10), internal::random<Scalar>(-1,2));
+      data2[i] = internal::random<Scalar>(0.01,1) * std::pow(
+          Scalar(10), internal::random<Scalar>(-1,2));
+  }
+
+#if EIGEN_HAS_C99_MATH && (__cplusplus > 199711L)
+  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
+  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
+  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc);
+#endif
+
+}
+
+namespace Eigen {
+namespace test {
+
+template<typename Scalar,typename PacketType, bool IsComplex, bool IsInteger>
+struct runall {
+  static void run() {
+    packetmath_real<Scalar,PacketType>();
+  }
+};
+
+}
+}
+
+EIGEN_DECLARE_TEST(special_packetmath)
+{
+  g_first_pass = true;
+  for(int i = 0; i < g_repeat; i++) {
+
+    CALL_SUBTEST_1( test::runner<float>::run() );
+    CALL_SUBTEST_2( test::runner<double>::run() );
+    g_first_pass = false;
+  }
+}