diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
index 5bce194..c3270e1 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -42,6 +42,134 @@
   typedef numext::int16_t type;
 };
 
+/* polevl (modified for Eigen)
+ *
+ *      Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * Scalar x, y, coef[N+1];
+ *
+ * y = polevl<decltype(x), N>( x, coef);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates polynomial of degree N:
+ *
+ *                     2          N
+ * y  =  C  + C x + C x  +...+ C x
+ *        0    1     2          N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C  , ..., coef[N] = C  .
+ *            N                   0
+ *
+ *  The function p1evl() assumes that coef[N] = 1.0 and is
+ * omitted from the array.  Its calling arguments are
+ * otherwise the same as polevl().
+ *
+ *
+ * The Eigen implementation is templatized.  For best speed, store
+ * coef as a const array (constexpr), e.g.
+ *
+ * const double coef[] = {1.0, 2.0, 3.0, ...};
+ *
+ */
+template <typename Packet, int N>
+struct ppolevl {
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x,
+                                                          const typename unpacket_traits<Packet>::type coeff[]) {
+    EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE);
+    return pmadd(ppolevl<Packet, N - 1>::run(x, coeff), x, pset1<Packet>(coeff[N]));
+  }
+};
+
+template <typename Packet>
+struct ppolevl<Packet, 0> {
+  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x,
+                                                          const typename unpacket_traits<Packet>::type coeff[]) {
+    EIGEN_UNUSED_VARIABLE(x);
+    return pset1<Packet>(coeff[0]);
+  }
+};
+
+/* chbevl (modified for Eigen)
+ *
+ *     Evaluate Chebyshev series
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * Scalar x, y, coef[N], chebevl();
+ *
+ * y = chbevl( x, coef, N );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates the series
+ *
+ *        N-1
+ *         - '
+ *  y  =   >   coef[i] T (x/2)
+ *         -            i
+ *        i=0
+ *
+ * of Chebyshev polynomials Ti at argument x/2.
+ *
+ * Coefficients are stored in reverse order, i.e. the zero
+ * order term is last in the array.  Note N is the number of
+ * coefficients, not the order.
+ *
+ * If coefficients are for the interval a to b, x must
+ * have been transformed to x -> 2(2x - b - a)/(b-a) before
+ * entering the routine.  This maps x from (a, b) to (-1, 1),
+ * over which the Chebyshev polynomials are defined.
+ *
+ * If the coefficients are for the inverted interval, in
+ * which (a, b) is mapped to (1/b, 1/a), the transformation
+ * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
+ * this becomes x -> 4a/x - 1.
+ *
+ *
+ *
+ * SPEED:
+ *
+ * Taking advantage of the recurrence properties of the
+ * Chebyshev polynomials, the routine requires one more
+ * addition per loop than evaluating a nested polynomial of
+ * the same degree.
+ *
+ */
+
+template <typename Packet, int N>
+struct pchebevl {
+  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Packet run(Packet x,
+                                                          const typename unpacket_traits<Packet>::type coef[]) {
+    typedef typename unpacket_traits<Packet>::type Scalar;
+    Packet b0 = pset1<Packet>(coef[0]);
+    Packet b1 = pset1<Packet>(static_cast<Scalar>(0.f));
+    Packet b2;
+
+    for (int i = 1; i < N; i++) {
+      b2 = b1;
+      b1 = b0;
+      b0 = psub(pmadd(x, b1, pset1<Packet>(coef[i])), b2);
+    }
+
+    return pmul(pset1<Packet>(static_cast<Scalar>(0.5f)), psub(b0, b2));
+  }
+};
+
 template <typename Packet>
 EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Packet pfrexp_generic_get_biased_exponent(const Packet& a) {
   typedef typename unpacket_traits<Packet>::type Scalar;
@@ -193,21 +321,14 @@
   e = psub(e, pand(cst_1, mask));
   x = padd(x, tmp);
 
-  // Polynomial coefficients for rational (3,3) r(x) = p(x)/q(x)
+  // Polynomial coefficients for rational r(x) = p(x)/q(x)
   // approximating log(1+x) on [sqrt(0.5)-1;sqrt(2)-1].
-  const Packet cst_p1 = pset1<Packet>(1.0000000190281136f);
-  const Packet cst_p2 = pset1<Packet>(1.0000000190281063f);
-  const Packet cst_p3 = pset1<Packet>(0.18256296349849254f);
-  const Packet cst_q1 = pset1<Packet>(1.4999999999999927f);
-  const Packet cst_q2 = pset1<Packet>(0.59923249590823520f);
-  const Packet cst_q3 = pset1<Packet>(0.049616247954120038f);
+  constexpr float alpha[] = {0.18256296349849254f, 1.0000000190281063f, 1.0000000190281136f};
+  constexpr float beta[] = {0.049616247954120038f, 0.59923249590823520f, 1.4999999999999927f, 1.0f};
 
-  Packet p = pmadd(x, cst_p3, cst_p2);
-  p = pmadd(x, p, cst_p1);
+  Packet p = ppolevl<Packet, 2>::run(x, alpha);
   p = pmul(x, p);
-  Packet q = pmadd(x, cst_q3, cst_q2);
-  q = pmadd(x, q, cst_q1);
-  q = pmadd(x, q, cst_1);
+  Packet q = ppolevl<Packet, 3>::run(x, beta);
   x = pdiv(p, q);
 
   // Add the logarithm of the exponent back to the result of the interpolation.
@@ -919,13 +1040,6 @@
   const Packet cst_one = pset1<Packet>(1.0f);
   const Packet cst_two = pset1<Packet>(2.0f);
   const Packet cst_pi_over_two = pset1<Packet>(kPiOverTwo);
-  // For |x| < 0.5 approximate asin(x)/x by an 8th order polynomial with
-  // even terms only.
-  const Packet p9 = pset1<Packet>(5.08838854730129241943359375e-2f);
-  const Packet p7 = pset1<Packet>(3.95139865577220916748046875e-2f);
-  const Packet p5 = pset1<Packet>(7.550220191478729248046875e-2f);
-  const Packet p3 = pset1<Packet>(0.16664917767047882080078125f);
-  const Packet p1 = pset1<Packet>(1.00000011920928955078125f);
 
   const Packet abs_x = pabs(x_in);
   const Packet sign_mask = pandnot(x_in, abs_x);
@@ -941,13 +1055,11 @@
   const Packet x = pselect(large_mask, x_large, abs_x);
   const Packet x2 = pmul(x, x);
 
-  // Compute polynomial.
-  // x * (p1 + x^2*(p3 + x^2*(p5 + x^2*(p7 + x^2*p9))))
-
-  Packet p = pmadd(p9, x2, p7);
-  p = pmadd(p, x2, p5);
-  p = pmadd(p, x2, p3);
-  p = pmadd(p, x2, p1);
+  // For |x| < 0.5 approximate asin(x)/x by an 8th order polynomial with
+  // even terms only.
+  constexpr float alpha[] = {5.08838854730129241943359375e-2f, 3.95139865577220916748046875e-2f,
+                             7.550220191478729248046875e-2f, 0.16664917767047882080078125f, 1.00000011920928955078125f};
+  Packet p = ppolevl<Packet, 4>::run(x2, alpha);
   p = pmul(p, x);
 
   const Packet p_large = pnmadd(cst_two, p, cst_pi_over_two);
@@ -967,34 +1079,17 @@
 template <>
 template <typename Packet>
 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patan_reduced<double>::run(const Packet& x) {
-  const Packet p1 = pset1<Packet>(3.3004361289279920e-01);
-  const Packet p3 = pset1<Packet>(8.2704055405494614e-01);
-  const Packet p5 = pset1<Packet>(7.5365702534987022e-01);
-  const Packet p7 = pset1<Packet>(3.0409318473444424e-01);
-  const Packet p9 = pset1<Packet>(5.2574296781008604e-02);
-  const Packet p11 = pset1<Packet>(3.0917513112462781e-03);
-  const Packet p13 = pset1<Packet>(2.6667153866462208e-05);
-  const Packet q0 = p1;
-  const Packet q2 = pset1<Packet>(9.3705509168587852e-01);
-  const Packet q4 = pset1<Packet>(1.0);
-  const Packet q6 = pset1<Packet>(4.9716458728465573e-01);
-  const Packet q8 = pset1<Packet>(1.1548932646420353e-01);
-  const Packet q10 = pset1<Packet>(1.0899150928962708e-02);
-  const Packet q12 = pset1<Packet>(2.7311202462436667e-04);
+  constexpr double alpha[] = {2.6667153866462208e-05, 3.0917513112462781e-03, 5.2574296781008604e-02,
+                              3.0409318473444424e-01, 7.5365702534987022e-01, 8.2704055405494614e-01,
+                              3.3004361289279920e-01};
+
+  constexpr double beta[] = {
+      2.7311202462436667e-04, 1.0899150928962708e-02, 1.1548932646420353e-01, 4.9716458728465573e-01, 1.0,
+      9.3705509168587852e-01, 3.3004361289279920e-01};
 
   Packet x2 = pmul(x, x);
-  Packet p = pmadd(p13, x2, p11);
-  Packet q = pmadd(q12, x2, q10);
-  p = pmadd(p, x2, p9);
-  q = pmadd(q, x2, q8);
-  p = pmadd(p, x2, p7);
-  q = pmadd(q, x2, q6);
-  p = pmadd(p, x2, p5);
-  q = pmadd(q, x2, q4);
-  p = pmadd(p, x2, p3);
-  q = pmadd(q, x2, q2);
-  p = pmadd(p, x2, p1);
-  q = pmadd(q, x2, q0);
+  Packet p = ppolevl<Packet, 6>::run(x2, alpha);
+  Packet q = ppolevl<Packet, 6>::run(x2, beta);
   return pmul(x, pdiv(p, q));
 }
 
@@ -1002,20 +1097,14 @@
 template <>
 template <typename Packet>
 EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patan_reduced<float>::run(const Packet& x) {
-  const Packet p1 = pset1<Packet>(8.109951019287109375e-01f);
-  const Packet p3 = pset1<Packet>(7.296695709228515625e-01f);
-  const Packet p5 = pset1<Packet>(1.12026982009410858154296875e-01f);
-  const Packet q0 = p1;
-  const Packet q2 = pset1<Packet>(1.0f);
-  const Packet q4 = pset1<Packet>(2.8318560123443603515625e-01f);
-  const Packet q6 = pset1<Packet>(1.00917108356952667236328125e-02f);
+  constexpr float alpha[] = {1.12026982009410858154296875e-01f, 7.296695709228515625e-01f, 8.109951019287109375e-01f};
+
+  constexpr float beta[] = {1.00917108356952667236328125e-02f, 2.8318560123443603515625e-01f, 1.0f,
+                            8.109951019287109375e-01f};
 
   Packet x2 = pmul(x, x);
-  Packet p = pmadd(p5, x2, p3);
-  Packet q = pmadd(q6, x2, q4);
-  p = pmadd(p, x2, p1);
-  q = pmadd(q, x2, q2);
-  q = pmadd(q, x2, q0);
+  Packet p = ppolevl<Packet, 2>::run(x2, alpha);
+  Packet q = ppolevl<Packet, 3>::run(x2, beta);
   return pmul(x, pdiv(p, q));
 }
 
@@ -1076,34 +1165,20 @@
   //   --output=tanhf.sollya --dispCoeff="dec"
 
   // The monomial coefficients of the numerator polynomial (odd).
-  const T alpha_3 = pset1<T>(1.340216100e-1f);
-  const T alpha_5 = pset1<T>(3.520756727e-3f);
-  const T alpha_7 = pset1<T>(2.102733560e-5f);
-  const T alpha_9 = pset1<T>(1.394553628e-8f);
+  constexpr float alpha[] = {1.394553628e-8f, 2.102733560e-5f, 3.520756727e-3f, 1.340216100e-1f};
 
   // The monomial coefficients of the denominator polynomial (even).
-  const T beta_0 = pset1<T>(1.0f);
-  const T beta_2 = pset1<T>(4.673548340e-1f);
-  const T beta_4 = pset1<T>(2.597254514e-2f);
-  const T beta_6 = pset1<T>(3.326951409e-4f);
-  const T beta_8 = pset1<T>(8.015776984e-7f);
+  constexpr float beta[] = {8.015776984e-7f, 3.326951409e-4f, 2.597254514e-2f, 4.673548340e-1f, 1.0f};
 
   // Since the polynomials are odd/even, we need x^2.
   const T x2 = pmul(x, x);
   const T x3 = pmul(x2, x);
 
-  // Interleave the evaluation of the numerator polynomial p and
-  // denominator polynomial q.
-  T p = pmadd(x2, alpha_9, alpha_7);
-  T q = pmadd(x2, beta_8, beta_6);
-  p = pmadd(x2, p, alpha_5);
-  q = pmadd(x2, q, beta_4);
-  p = pmadd(x2, p, alpha_3);
-  q = pmadd(x2, q, beta_2);
-  // Take advantage of the fact that alpha_1 = 1 to compute
-  // x*(x^2*p + alpha_1) = x^3 * p + x.
+  T p = ppolevl<T, 3>::run(x2, alpha);
+  T q = ppolevl<T, 4>::run(x2, beta);
+  // Take advantage of the fact that the constant term in p is 1 to compute
+  // x*(x^2*p + 1) = x^3 * p + x.
   p = pmadd(x3, p, x);
-  q = pmadd(x2, q, beta_0);
 
   // Divide the numerator by the denominator.
   return pdiv(p, q);
@@ -1141,27 +1216,16 @@
   //   --denF="[D]" --log --output=tanh.sollya --dispCoeff="dec"
 
   // The monomial coefficients of the numerator polynomial (odd).
-  const T alpha_3 = pset1<T>(1.5184719640284322e-01);
-  const T alpha_5 = pset1<T>(5.9809711724441161e-03);
-  const T alpha_7 = pset1<T>(9.3839087674268880e-05);
-  const T alpha_9 = pset1<T>(6.8644367682497074e-07);
-  const T alpha_11 = pset1<T>(2.4618379131293676e-09);
-  const T alpha_13 = pset1<T>(4.2303918148209176e-12);
-  const T alpha_15 = pset1<T>(3.1309488231386680e-15);
-  const T alpha_17 = pset1<T>(7.6534862268749319e-19);
-  const T alpha_19 = pset1<T>(2.6158007860482230e-23);
+  constexpr double alpha[] = {2.6158007860482230e-23, 7.6534862268749319e-19, 3.1309488231386680e-15,
+                              4.2303918148209176e-12, 2.4618379131293676e-09, 6.8644367682497074e-07,
+                              9.3839087674268880e-05, 5.9809711724441161e-03, 1.5184719640284322e-01};
 
   // The monomial coefficients of the denominator polynomial (even).
-  const T beta_0 = pset1<T>(1.0);
-  const T beta_2 = pset1<T>(4.851805297361760360e-01);
-  const T beta_4 = pset1<T>(3.437448108450402717e-02);
-  const T beta_6 = pset1<T>(8.295161192716231542e-04);
-  const T beta_8 = pset1<T>(8.785185266237658698e-06);
-  const T beta_10 = pset1<T>(4.492975677839633985e-08);
-  const T beta_12 = pset1<T>(1.123643448069621992e-10);
-  const T beta_14 = pset1<T>(1.293019623712687916e-13);
-  const T beta_16 = pset1<T>(5.782506856739003571e-17);
-  const T beta_18 = pset1<T>(6.463747022670968018e-21);
+  constexpr double beta[] = {6.463747022670968018e-21, 5.782506856739003571e-17,
+                             1.293019623712687916e-13, 1.123643448069621992e-10,
+                             4.492975677839633985e-08, 8.785185266237658698e-06,
+                             8.295161192716231542e-04, 3.437448108450402717e-02,
+                             4.851805297361760360e-01, 1.0};
 
   // Since the polynomials are odd/even, we need x^2.
   const T x2 = pmul(x, x);
@@ -1169,26 +1233,11 @@
 
   // Interleave the evaluation of the numerator polynomial p and
   // denominator polynomial q.
-  T p = pmadd(x2, alpha_19, alpha_17);
-  T q = pmadd(x2, beta_18, beta_16);
-  p = pmadd(x2, p, alpha_15);
-  q = pmadd(x2, q, beta_14);
-  p = pmadd(x2, p, alpha_13);
-  q = pmadd(x2, q, beta_12);
-  p = pmadd(x2, p, alpha_11);
-  q = pmadd(x2, q, beta_10);
-  p = pmadd(x2, p, alpha_9);
-  q = pmadd(x2, q, beta_8);
-  p = pmadd(x2, p, alpha_7);
-  q = pmadd(x2, q, beta_6);
-  p = pmadd(x2, p, alpha_5);
-  q = pmadd(x2, q, beta_4);
-  p = pmadd(x2, p, alpha_3);
-  q = pmadd(x2, q, beta_2);
-  // Take advantage of the fact that alpha_1 = 1 to compute
-  // x*(x^2*p + alpha_1) = x^3 * p + x.
+  T p = ppolevl<T, 8>::run(x2, alpha);
+  T q = ppolevl<T, 9>::run(x2, beta);
+  // Take advantage of the fact that the constant term in p is 1 to compute
+  // x*(x^2*p + 1) = x^3 * p + x.
   p = pmadd(x3, p, x);
-  q = pmadd(x2, q, beta_0);
 
   // Divide the numerator by the denominator.
   return pdiv(p, q);
@@ -1200,18 +1249,13 @@
   static_assert(std::is_same<Scalar, float>::value, "Scalar type must be float");
 
   // For |x| in [0:0.5] we use a polynomial approximation of the form
-  // P(x) = x + x^3*(c3 + x^2 * (c5 + x^2 * (... x^2 * c11) ... )).
-  const Packet C3 = pset1<Packet>(0.3333373963832855224609375f);
-  const Packet C5 = pset1<Packet>(0.1997792422771453857421875f);
-  const Packet C7 = pset1<Packet>(0.14672131836414337158203125f);
-  const Packet C9 = pset1<Packet>(8.2311116158962249755859375e-2f);
-  const Packet C11 = pset1<Packet>(0.1819281280040740966796875f);
+  // P(x) = x + x^3*(alpha[4] + x^2 * (alpha[3] + x^2 * (... x^2 * alpha[0]) ... )).
+  constexpr float alpha[] = {0.1819281280040740966796875f, 8.2311116158962249755859375e-2f,
+                             0.14672131836414337158203125f, 0.1997792422771453857421875f, 0.3333373963832855224609375f};
   const Packet x2 = pmul(x, x);
-  Packet p = pmadd(C11, x2, C9);
-  p = pmadd(x2, p, C7);
-  p = pmadd(x2, p, C5);
-  p = pmadd(x2, p, C3);
-  p = pmadd(pmul(x, x2), p, x);
+  const Packet x3 = pmul(x, x2);
+  Packet p = ppolevl<Packet, 4>::run(x2, alpha);
+  p = pmadd(x3, p, x);
 
   // For |x| in ]0.5:1.0] we use atanh = 0.5*ln((1+x)/(1-x));
   const Packet half = pset1<Packet>(0.5f);
@@ -1234,30 +1278,16 @@
   static_assert(std::is_same<Scalar, double>::value, "Scalar type must be double");
   // For x in [-0.5:0.5] we use a rational approximation of the form
   // R(x) = x + x^3*P(x^2)/Q(x^2), where P is or order 4 and Q is of order 5.
-  const Packet p0 = pset1<Packet>(1.2306328729812676e-01);
-  const Packet p2 = pset1<Packet>(-2.5949536095445679e-01);
-  const Packet p4 = pset1<Packet>(1.8185306179826699e-01);
-  const Packet p6 = pset1<Packet>(-4.7129526768798737e-02);
-  const Packet p8 = pset1<Packet>(3.3071338469301391e-03);
+  constexpr double alpha[] = {3.3071338469301391e-03, -4.7129526768798737e-02, 1.8185306179826699e-01,
+                              -2.5949536095445679e-01, 1.2306328729812676e-01};
 
-  const Packet q0 = pset1<Packet>(3.6918986189438030e-01);
-  const Packet q2 = pset1<Packet>(-1.0000000000000000e+00);
-  const Packet q4 = pset1<Packet>(9.8733495886883648e-01);
-  const Packet q6 = pset1<Packet>(-4.2828141436397615e-01);
-  const Packet q8 = pset1<Packet>(7.6391885763341910e-02);
-  const Packet q10 = pset1<Packet>(-3.8679974580640881e-03);
+  constexpr double beta[] = {-3.8679974580640881e-03, 7.6391885763341910e-02,  -4.2828141436397615e-01,
+                             9.8733495886883648e-01,  -1.0000000000000000e+00, 3.6918986189438030e-01};
 
   const Packet x2 = pmul(x, x);
   const Packet x3 = pmul(x, x2);
-  Packet q = pmadd(q10, x2, q8);
-  Packet p = pmadd(p8, x2, p6);
-  q = pmadd(x2, q, q6);
-  p = pmadd(x2, p, p4);
-  q = pmadd(x2, q, q4);
-  p = pmadd(x2, p, p2);
-  q = pmadd(x2, q, q2);
-  p = pmadd(x2, p, p0);
-  q = pmadd(x2, q, q0);
+  Packet p = ppolevl<Packet, 4>::run(x2, alpha);
+  Packet q = ppolevl<Packet, 5>::run(x2, beta);
   Packet y_small = pmadd(x3, pdiv(p, q), x);
 
   // For |x| in ]0.5:1.0] we use atanh = 0.5*ln((1+x)/(1-x));
@@ -2200,134 +2230,6 @@
                                                          pselect(negate_pow_abs, pnegate(pow_abs), pow_abs)))))));
 }
 
-/* polevl (modified for Eigen)
- *
- *      Evaluate polynomial
- *
- *
- *
- * SYNOPSIS:
- *
- * int N;
- * Scalar x, y, coef[N+1];
- *
- * y = polevl<decltype(x), N>( x, coef);
- *
- *
- *
- * DESCRIPTION:
- *
- * Evaluates polynomial of degree N:
- *
- *                     2          N
- * y  =  C  + C x + C x  +...+ C x
- *        0    1     2          N
- *
- * Coefficients are stored in reverse order:
- *
- * coef[0] = C  , ..., coef[N] = C  .
- *            N                   0
- *
- *  The function p1evl() assumes that coef[N] = 1.0 and is
- * omitted from the array.  Its calling arguments are
- * otherwise the same as polevl().
- *
- *
- * The Eigen implementation is templatized.  For best speed, store
- * coef as a const array (constexpr), e.g.
- *
- * const double coef[] = {1.0, 2.0, 3.0, ...};
- *
- */
-template <typename Packet, int N>
-struct ppolevl {
-  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x,
-                                                          const typename unpacket_traits<Packet>::type coeff[]) {
-    EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE);
-    return pmadd(ppolevl<Packet, N - 1>::run(x, coeff), x, pset1<Packet>(coeff[N]));
-  }
-};
-
-template <typename Packet>
-struct ppolevl<Packet, 0> {
-  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& x,
-                                                          const typename unpacket_traits<Packet>::type coeff[]) {
-    EIGEN_UNUSED_VARIABLE(x);
-    return pset1<Packet>(coeff[0]);
-  }
-};
-
-/* chbevl (modified for Eigen)
- *
- *     Evaluate Chebyshev series
- *
- *
- *
- * SYNOPSIS:
- *
- * int N;
- * Scalar x, y, coef[N], chebevl();
- *
- * y = chbevl( x, coef, N );
- *
- *
- *
- * DESCRIPTION:
- *
- * Evaluates the series
- *
- *        N-1
- *         - '
- *  y  =   >   coef[i] T (x/2)
- *         -            i
- *        i=0
- *
- * of Chebyshev polynomials Ti at argument x/2.
- *
- * Coefficients are stored in reverse order, i.e. the zero
- * order term is last in the array.  Note N is the number of
- * coefficients, not the order.
- *
- * If coefficients are for the interval a to b, x must
- * have been transformed to x -> 2(2x - b - a)/(b-a) before
- * entering the routine.  This maps x from (a, b) to (-1, 1),
- * over which the Chebyshev polynomials are defined.
- *
- * If the coefficients are for the inverted interval, in
- * which (a, b) is mapped to (1/b, 1/a), the transformation
- * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
- * this becomes x -> 4a/x - 1.
- *
- *
- *
- * SPEED:
- *
- * Taking advantage of the recurrence properties of the
- * Chebyshev polynomials, the routine requires one more
- * addition per loop than evaluating a nested polynomial of
- * the same degree.
- *
- */
-
-template <typename Packet, int N>
-struct pchebevl {
-  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Packet run(Packet x,
-                                                          const typename unpacket_traits<Packet>::type coef[]) {
-    typedef typename unpacket_traits<Packet>::type Scalar;
-    Packet b0 = pset1<Packet>(coef[0]);
-    Packet b1 = pset1<Packet>(static_cast<Scalar>(0.f));
-    Packet b2;
-
-    for (int i = 1; i < N; i++) {
-      b2 = b1;
-      b1 = b0;
-      b0 = psub(pmadd(x, b1, pset1<Packet>(coef[i])), b2);
-    }
-
-    return pmul(pset1<Packet>(static_cast<Scalar>(0.5f)), psub(b0, b2));
-  }
-};
-
 namespace unary_pow {
 
 template <typename ScalarExponent, bool IsInteger = NumTraits<ScalarExponent>::IsInteger>
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
index 56615d3..9f12a34 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
@@ -131,8 +131,8 @@
 template <>
 struct digamma_impl_maybe_poly<float> {
   EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE float run(const float s) {
-    const float A[] = {-4.16666666666666666667E-3f, 3.96825396825396825397E-3f, -8.33333333333333333333E-3f,
-                       8.33333333333333333333E-2f};
+    constexpr float A[] = {-4.16666666666666666667E-3f, 3.96825396825396825397E-3f, -8.33333333333333333333E-3f,
+                           8.33333333333333333333E-2f};
 
     float z;
     if (s < 1.0e8f) {
@@ -146,9 +146,9 @@
 template <>
 struct digamma_impl_maybe_poly<double> {
   EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE double run(const double s) {
-    const double A[] = {8.33333333333333333333E-2,  -2.10927960927960927961E-2, 7.57575757575757575758E-3,
-                        -4.16666666666666666667E-3, 3.96825396825396825397E-3,  -8.33333333333333333333E-3,
-                        8.33333333333333333333E-2};
+    constexpr double A[] = {8.33333333333333333333E-2,  -2.10927960927960927961E-2, 7.57575757575757575758E-3,
+                            -4.16666666666666666667E-3, 3.96825396825396825397E-3,  -8.33333333333333333333E-3,
+                            8.33333333333333333333E-2};
 
     double z;
     if (s < 1.0e17) {
@@ -282,20 +282,14 @@
 template <typename T>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erf_float(const T& x) {
   // The monomial coefficients of the numerator polynomial (odd).
-  const T alpha_1 = pset1<T>(1.128379106521606445312500000000000000e+00f);
-  const T alpha_3 = pset1<T>(1.874160766601562500000000000000000000e-01f);
-  const T alpha_5 = pset1<T>(5.243302136659622192382812500000000000e-02f);
-  const T alpha_7 = pset1<T>(3.658048342913389205932617187500000000e-03f);
-  const T alpha_9 = pset1<T>(2.861979592125862836837768554687500000e-04f);
-  const T alpha_11 = pset1<T>(2.123732201653183437883853912353515625e-06f);
+  constexpr float alpha[] = {2.123732201653183437883853912353515625e-06f, 2.861979592125862836837768554687500000e-04f,
+                             3.658048342913389205932617187500000000e-03f, 5.243302136659622192382812500000000000e-02f,
+                             1.874160766601562500000000000000000000e-01f, 1.128379106521606445312500000000000000e+00f};
 
   // The monomial coefficients of the denominator polynomial (even).
-  const T beta_0 = pset1<T>(1.0f);
-  const T beta_2 = pset1<T>(4.99425798654556274414062500000e-01f);
-  const T beta_4 = pset1<T>(1.12945675849914550781250000000e-01f);
-  const T beta_6 = pset1<T>(1.47520881146192550659179687500e-02f);
-  const T beta_8 = pset1<T>(1.14329601638019084930419921875e-03f);
-  const T beta_10 = pset1<T>(3.89185734093189239501953125000e-05f);
+  constexpr float beta[] = {3.89185734093189239501953125000e-05f, 1.14329601638019084930419921875e-03f,
+                            1.47520881146192550659179687500e-02f, 1.12945675849914550781250000000e-01f,
+                            4.99425798654556274414062500000e-01f, 1.0f};
 
   // Since the polynomials are odd/even, we need x^2.
   // Since erf(4) == 1 in float, we clamp x^2 to 16 to avoid
@@ -303,19 +297,11 @@
   const T x2 = pmin(pset1<T>(16.0f), pmul(x, x));
 
   // Evaluate the numerator polynomial p.
-  T p = pmadd(x2, alpha_11, alpha_9);
-  p = pmadd(x2, p, alpha_7);
-  p = pmadd(x2, p, alpha_5);
-  p = pmadd(x2, p, alpha_3);
-  p = pmadd(x2, p, alpha_1);
+  T p = ppolevl<T, 5>::run(x2, alpha);
   p = pmul(x, p);
 
   // Evaluate the denominator polynomial p.
-  T q = pmadd(x2, beta_10, beta_8);
-  q = pmadd(x2, q, beta_6);
-  q = pmadd(x2, q, beta_4);
-  q = pmadd(x2, q, beta_2);
-  q = pmadd(x2, q, beta_0);
+  T q = ppolevl<T, 5>::run(x2, beta);
   const T r = pdiv(p, q);
 
   // Clamp to [-1:1].
@@ -475,18 +461,18 @@
 
 template <typename T, typename ScalarType>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_ndtri_gt_exp_neg_two(const T& b) {
-  const ScalarType p0[] = {ScalarType(-5.99633501014107895267e1), ScalarType(9.80010754185999661536e1),
-                           ScalarType(-5.66762857469070293439e1), ScalarType(1.39312609387279679503e1),
-                           ScalarType(-1.23916583867381258016e0)};
-  const ScalarType q0[] = {ScalarType(1.0),
-                           ScalarType(1.95448858338141759834e0),
-                           ScalarType(4.67627912898881538453e0),
-                           ScalarType(8.63602421390890590575e1),
-                           ScalarType(-2.25462687854119370527e2),
-                           ScalarType(2.00260212380060660359e2),
-                           ScalarType(-8.20372256168333339912e1),
-                           ScalarType(1.59056225126211695515e1),
-                           ScalarType(-1.18331621121330003142e0)};
+  constexpr ScalarType p0[] = {ScalarType(-5.99633501014107895267e1), ScalarType(9.80010754185999661536e1),
+                               ScalarType(-5.66762857469070293439e1), ScalarType(1.39312609387279679503e1),
+                               ScalarType(-1.23916583867381258016e0)};
+  constexpr ScalarType q0[] = {ScalarType(1.0),
+                               ScalarType(1.95448858338141759834e0),
+                               ScalarType(4.67627912898881538453e0),
+                               ScalarType(8.63602421390890590575e1),
+                               ScalarType(-2.25462687854119370527e2),
+                               ScalarType(2.00260212380060660359e2),
+                               ScalarType(-8.20372256168333339912e1),
+                               ScalarType(1.59056225126211695515e1),
+                               ScalarType(-1.18331621121330003142e0)};
   const T sqrt2pi = pset1<T>(ScalarType(2.50662827463100050242e0));
   const T half = pset1<T>(ScalarType(0.5));
   T c, c2, ndtri_gt_exp_neg_two;
@@ -503,20 +489,20 @@
   /* Approximation for interval z = sqrt(-2 log a ) between 2 and 8
    * i.e., a between exp(-2) = .135 and exp(-32) = 1.27e-14.
    */
-  const ScalarType p1[] = {ScalarType(4.05544892305962419923e0),   ScalarType(3.15251094599893866154e1),
-                           ScalarType(5.71628192246421288162e1),   ScalarType(4.40805073893200834700e1),
-                           ScalarType(1.46849561928858024014e1),   ScalarType(2.18663306850790267539e0),
-                           ScalarType(-1.40256079171354495875e-1), ScalarType(-3.50424626827848203418e-2),
-                           ScalarType(-8.57456785154685413611e-4)};
-  const ScalarType q1[] = {ScalarType(1.0),
-                           ScalarType(1.57799883256466749731e1),
-                           ScalarType(4.53907635128879210584e1),
-                           ScalarType(4.13172038254672030440e1),
-                           ScalarType(1.50425385692907503408e1),
-                           ScalarType(2.50464946208309415979e0),
-                           ScalarType(-1.42182922854787788574e-1),
-                           ScalarType(-3.80806407691578277194e-2),
-                           ScalarType(-9.33259480895457427372e-4)};
+  constexpr ScalarType p1[] = {ScalarType(4.05544892305962419923e0),   ScalarType(3.15251094599893866154e1),
+                               ScalarType(5.71628192246421288162e1),   ScalarType(4.40805073893200834700e1),
+                               ScalarType(1.46849561928858024014e1),   ScalarType(2.18663306850790267539e0),
+                               ScalarType(-1.40256079171354495875e-1), ScalarType(-3.50424626827848203418e-2),
+                               ScalarType(-8.57456785154685413611e-4)};
+  constexpr ScalarType q1[] = {ScalarType(1.0),
+                               ScalarType(1.57799883256466749731e1),
+                               ScalarType(4.53907635128879210584e1),
+                               ScalarType(4.13172038254672030440e1),
+                               ScalarType(1.50425385692907503408e1),
+                               ScalarType(2.50464946208309415979e0),
+                               ScalarType(-1.42182922854787788574e-1),
+                               ScalarType(-3.80806407691578277194e-2),
+                               ScalarType(-9.33259480895457427372e-4)};
   /* Approximation for interval z = sqrt(-2 log a ) between 8 and 64
    * i.e., a between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    */
@@ -1267,7 +1253,7 @@
     int i;
     Scalar p, r, a, b, k, s, t, w;
 
-    const Scalar A[] = {
+    constexpr Scalar A[] = {
         Scalar(12.0),
         Scalar(-720.0),
         Scalar(30240.0),