--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h	
+++ a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h	
@@ -228,16 +228,20 @@ Packet pexp_float(const Packet _x)
   y = pmadd(y, r, cst_cephes_exp_p5);
   y = pmadd(y, r2, r);
   y = padd(y, cst_1);
 
   // Return 2^m * exp(r).
   return pmax(pldexp(y,m), _x);
 }
 
+// make it the default path for scalar float
+template<>
+float pexp(const float& a) { return pexp_float(a); }
+
 template <typename Packet>
 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 EIGEN_UNUSED
 Packet pexp_double(const Packet _x)
 {
   Packet x = _x;
 
   const Packet cst_1 = pset1<Packet>(1.0);
@@ -296,16 +300,20 @@ Packet pexp_double(const Packet _x)
   x = pdiv(px, psub(qx, px));
   x = pmadd(cst_2, x, cst_1);
 
   // Construct the result 2^n * exp(g) = e * x. The max is used to catch
   // non-finite values in the input.
   return pmax(pldexp(x,fx), _x);
 }
 
+// make it the default path for scalar double
+template<>
+double pexp(const double& a) { return pexp_double(a); }
+
 // The following code is inspired by the following stack-overflow answer:
 //   https://stackoverflow.com/questions/30463616/payne-hanek-algorithm-implementation-in-c/30465751#30465751
 // It has been largely optimized:
 //  - By-pass calls to frexp.
 //  - Aligned loads of required 96 bits of 2/pi. This is accomplished by
 //    (1) balancing the mantissa and exponent to the required bits of 2/pi are
 //    aligned on 8-bits, and (2) replicating the storage of the bits of 2/pi.
 //  - Avoid a branch in rounding and extraction of the remaining fractional part.
--- a/Eigen/src/Core/functors/UnaryFunctors.h	
+++ a/Eigen/src/Core/functors/UnaryFunctors.h	
@@ -889,18 +889,17 @@ struct functor_traits<scalar_sign_op<Sca
 /** \internal
   * \brief Template functor to compute the logistic function of a scalar
   * \sa class CwiseUnaryOp, ArrayBase::logistic()
   */
 template <typename T>
 struct scalar_logistic_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(const T& x) const {
-    const T one = T(1);
-    return one / (one + numext::exp(-x));
+    return packetOp(x);
   }
 
   template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
   Packet packetOp(const Packet& x) const {
     const Packet one = pset1<Packet>(T(1));
     return pdiv(one, padd(one, pexp(pnegate(x))));
   }
 };
@@ -914,19 +913,17 @@ struct scalar_logistic_op {
   *  logistic is interpolated because it was easier to make the fit converge.
   *
   */
 
 template <>
 struct scalar_logistic_op<float> {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float operator()(const float& x) const {
-    if (x < -18.0f) return 0.0f;
-    else if (x > 18.0f) return 1.0f;
-    else return 1.0f / (1.0f + numext::exp(-x));
+    return packetOp(x);
   }
 
   template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
   Packet packetOp(const Packet& _x) const {
     // Clamp the inputs to the range [-18, 18] since anything outside
     // this range is 0.0f or 1.0f in single-precision.
     const Packet x = pmax(pmin(_x, pset1<Packet>(18.0)), pset1<Packet>(-18.0));
 
--- a/test/packetmath.cpp	
+++ a/test/packetmath.cpp	
@@ -566,16 +566,31 @@ template<typename Scalar,typename Packet
   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;
     h.store(data2, internal::ptanh(h.load(data1)));
     VERIFY((numext::isnan)(data2[0]));
   }
 
+  {
+    internal::scalar_logistic_op<Scalar> logistic;
+    for (int i=0; i<size; ++i)
+    {
+      data1[i] = internal::random<Scalar>(-20,20);
+    }
+    internal::pstore(data2, logistic.packetOp(internal::pload<Packet>(data1)));
+    for (int i=0; i<PacketSize; ++i) {
+      VERIFY_IS_APPROX(data2[i],logistic(data1[i]));
+      #ifdef EIGEN_VECTORIZE // don't check for exactness when using the i387 FPU
+      VERIFY_IS_EQUAL(data2[i],logistic(data1[i]));
+      #endif
+    }
+  }
+
 #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) {
@@ -960,17 +975,16 @@ struct runner<Scalar,PacketType,false,fa
     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_6(( packetmath<half,internal::packet_traits<half>::type>() ));
     g_first_pass = false;
   }