bug #1232: refactor special functions as a new SpecialFunctions module, currently in unsupported/.

10 files changed, 102 insertions, 2011 deletions

diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h
index 76a75dee1..16a122370 100644
--- a/Eigen/src/Core/GenericPacketMath.h
+++ b/Eigen/src/Core/GenericPacketMath.h
@@ -435,42 +435,6 @@ Packet pfloor(const Packet& a) { using numext::floor; return floor(a); }
 template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 Packet pceil(const Packet& a) { using numext::ceil; return ceil(a); }
 
-/** \internal \returns the ln(|gamma(\a a)|) (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet plgamma(const Packet& a) { using numext::lgamma; return lgamma(a); }
-
-/** \internal \returns the derivative of lgamma, psi(\a a) (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet pdigamma(const Packet& a) { using numext::digamma; return digamma(a); }
-    
-/** \internal \returns the zeta function of two arguments (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet pzeta(const Packet& x, const Packet& q) { using numext::zeta; return zeta(x, q); }
-
-/** \internal \returns the polygamma function (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet ppolygamma(const Packet& n, const Packet& x) { using numext::polygamma; return polygamma(n, x); }
-
-/** \internal \returns the erf(\a a) (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet perf(const Packet& a) { using numext::erf; return erf(a); }
-
-/** \internal \returns the erfc(\a a) (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet perfc(const Packet& a) { using numext::erfc; return erfc(a); }
-
-/** \internal \returns the incomplete gamma function igamma(\a a, \a x) */
-template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
-Packet pigamma(const Packet& a, const Packet& x) { using numext::igamma; return igamma(a, x); }
-
-/** \internal \returns the complementary incomplete gamma function igammac(\a a, \a x) */
-template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
-Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; return igammac(a, x); }
-
-/** \internal \returns the complementary incomplete gamma function betainc(\a a, \a b, \a x) */
-template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
-Packet pbetainc(const Packet& a, const Packet& b,const Packet& x) { using numext::betainc; return betainc(a, b, x); }
-
 /***************************************************************************
 * The following functions might not have to be overwritten for vectorized types
 ***************************************************************************/
diff --git a/Eigen/src/Core/GlobalFunctions.h b/Eigen/src/Core/GlobalFunctions.h
index b9c3ec25b..879a93e6b 100644
--- a/Eigen/src/Core/GlobalFunctions.h
+++ b/Eigen/src/Core/GlobalFunctions.h
@@ -174,111 +174,6 @@ namespace Eigen
   }
 #endif
 
-  /** \cpp11 \returns an expression of the coefficient-wise igamma(\a a, \a x) to the given arrays.
-    *
-    * This function computes the coefficient-wise incomplete gamma function.
-    *
-    * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-    * or float/double in non c++11 mode, the user has to provide implementations of igammac(T,T) for any scalar
-    * type T to be supported.
-    *
-    * \sa Eigen::igammac(), Eigen::lgamma()
-    */
-  template<typename Derived,typename ExponentDerived>
-  inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_op<typename Derived::Scalar>, const Derived, const ExponentDerived>
-  igamma(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) 
-  {
-    return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_op<typename Derived::Scalar>, const Derived, const ExponentDerived>(
-      a.derived(),
-      x.derived()
-    );
-  }
-
-  /** \cpp11 \returns an expression of the coefficient-wise igammac(\a a, \a x) to the given arrays.
-    *
-    * This function computes the coefficient-wise complementary incomplete gamma function.
-    *
-    * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-    * or float/double in non c++11 mode, the user has to provide implementations of igammac(T,T) for any scalar
-    * type T to be supported.
-    *
-    * \sa Eigen::igamma(), Eigen::lgamma()
-    */
-  template<typename Derived,typename ExponentDerived>
-  inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igammac_op<typename Derived::Scalar>, const Derived, const ExponentDerived>
-  igammac(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) 
-  {
-    return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igammac_op<typename Derived::Scalar>, const Derived, const ExponentDerived>(
-      a.derived(),
-      x.derived()
-    );
-  }
-
-  /** \cpp11 \returns an expression of the coefficient-wise polygamma(\a n, \a x) to the given arrays.
-    *
-    * It returns the \a n -th derivative of the digamma(psi) evaluated at \c x.
-    *
-    * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-    * or float/double in non c++11 mode, the user has to provide implementations of polygamma(T,T) for any scalar
-    * type T to be supported.
-    *
-    * \sa Eigen::digamma()
-    */
-  // * \warning Be careful with the order of the parameters: x.polygamma(n) is equivalent to polygamma(n,x)
-  // * \sa ArrayBase::polygamma()
-  template<typename DerivedN,typename DerivedX>
-  inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_polygamma_op<typename DerivedX::Scalar>, const DerivedN, const DerivedX>
-  polygamma(const Eigen::ArrayBase<DerivedN>& n, const Eigen::ArrayBase<DerivedX>& x)
-  {
-    return Eigen::CwiseBinaryOp<Eigen::internal::scalar_polygamma_op<typename DerivedX::Scalar>, const DerivedN, const DerivedX>(
-      n.derived(),
-      x.derived()
-    );
-  }
-
-  /** \cpp11 \returns an expression of the coefficient-wise betainc(\a x, \a a, \a b) to the given arrays.
-    *
-    * This function computes the regularized incomplete beta function (integral).
-    *
-    * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-    * or float/double in non c++11 mode, the user has to provide implementations of betainc(T,T,T) for any scalar
-    * type T to be supported.
-    *
-    * \sa Eigen::betainc(), Eigen::lgamma()
-    */
-  template<typename ArgADerived, typename ArgBDerived, typename ArgXDerived>
-  inline const Eigen::CwiseTernaryOp<Eigen::internal::scalar_betainc_op<typename ArgXDerived::Scalar>, const ArgADerived, const ArgBDerived, const ArgXDerived>
-  betainc(const Eigen::ArrayBase<ArgADerived>& a, const Eigen::ArrayBase<ArgBDerived>& b, const Eigen::ArrayBase<ArgXDerived>& x)
-  {
-    return Eigen::CwiseTernaryOp<Eigen::internal::scalar_betainc_op<typename ArgXDerived::Scalar>, const ArgADerived, const ArgBDerived, const ArgXDerived>(
-      a.derived(),
-      b.derived(),
-      x.derived()
-    );
-  }
-
-
-  /** \returns an expression of the coefficient-wise zeta(\a x, \a q) to the given arrays.
-    *
-    * It returns the Riemann zeta function of two arguments \a x and \a q:
-    *
-    * \param x is the exposent, it must be > 1
-    * \param q is the shift, it must be > 0
-    *
-    * \note This function supports only float and double scalar types. To support other scalar types, the user has
-    * to provide implementations of zeta(T,T) for any scalar type T to be supported.
-    *
-    * \sa ArrayBase::zeta()
-    */
-  template<typename DerivedX,typename DerivedQ>
-  inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_zeta_op<typename DerivedX::Scalar>, const DerivedX, const DerivedQ>
-  zeta(const Eigen::ArrayBase<DerivedX>& x, const Eigen::ArrayBase<DerivedQ>& q)
-  {
-    return Eigen::CwiseBinaryOp<Eigen::internal::scalar_zeta_op<typename DerivedX::Scalar>, const DerivedX, const DerivedQ>(
-      x.derived(),
-      q.derived()
-    );
-  }
 
   namespace internal
   {
diff --git a/Eigen/src/Core/SpecialFunctions.h b/Eigen/src/Core/SpecialFunctions.h
deleted file mode 100644
index a657cb854..000000000
--- a/Eigen/src/Core/SpecialFunctions.h
+++ /dev/null
@@ -1,1551 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2015 Eugene Brevdo <ebrevdo@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/.
-
-#ifndef EIGEN_SPECIAL_FUNCTIONS_H
-#define EIGEN_SPECIAL_FUNCTIONS_H
-
-namespace Eigen {
-namespace internal {
-
-//  Parts of this code are based on the Cephes Math Library.
-//
-//  Cephes Math Library Release 2.8:  June, 2000
-//  Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
-//
-//  Permission has been kindly provided by the original author
-//  to incorporate the Cephes software into the Eigen codebase:
-//
-//    From: Stephen Moshier
-//    To: Eugene Brevdo
-//    Subject: Re: Permission to wrap several cephes functions in Eigen
-//
-//    Hello Eugene,
-//
-//    Thank you for writing.
-//
-//    If your licensing is similar to BSD, the formal way that has been
-//    handled is simply to add a statement to the effect that you are incorporating
-//    the Cephes software by permission of the author.
-//
-//    Good luck with your project,
-//    Steve
-
-namespace cephes {
-
-/* 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 Scalar, int N>
-struct polevl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar x, const Scalar coef[]) {
-    EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE);
-
-    return polevl<Scalar, N - 1>::run(x, coef) * x + coef[N];
-  }
-};
-
-template <typename Scalar>
-struct polevl<Scalar, 0> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar, const Scalar coef[]) {
-    return coef[0];
-  }
-};
-
-}  // end namespace cephes
-
-/****************************************************************************
- * Implementation of lgamma, requires C++11/C99                             *
- ****************************************************************************/
-
-template <typename Scalar>
-struct lgamma_impl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-template <typename Scalar>
-struct lgamma_retval {
-  typedef Scalar type;
-};
-
-#if EIGEN_HAS_C99_MATH
-template <>
-struct lgamma_impl<float> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float run(float x) { return ::lgammaf(x); }
-};
-
-template <>
-struct lgamma_impl<double> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double run(double x) { return ::lgamma(x); }
-};
-#endif
-
-/****************************************************************************
- * Implementation of digamma (psi), based on Cephes                         *
- ****************************************************************************/
-
-template <typename Scalar>
-struct digamma_retval {
-  typedef Scalar type;
-};
-
-/*
- *
- * Polynomial evaluation helper for the Psi (digamma) function.
- *
- * digamma_impl_maybe_poly::run(s) evaluates the asymptotic Psi expansion for
- * input Scalar s, assuming s is above 10.0.
- *
- * If s is above a certain threshold for the given Scalar type, zero
- * is returned.  Otherwise the polynomial is evaluated with enough
- * coefficients for results matching Scalar machine precision.
- *
- *
- */
-template <typename Scalar>
-struct digamma_impl_maybe_poly {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-
-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
-    };
-
-    float z;
-    if (s < 1.0e8f) {
-      z = 1.0f / (s * s);
-      return z * cephes::polevl<float, 3>::run(z, A);
-    } else return 0.0f;
-  }
-};
-
-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
-    };
-
-    double z;
-    if (s < 1.0e17) {
-      z = 1.0 / (s * s);
-      return z * cephes::polevl<double, 6>::run(z, A);
-    }
-    else return 0.0;
-  }
-};
-
-template <typename Scalar>
-struct digamma_impl {
-  EIGEN_DEVICE_FUNC
-  static Scalar run(Scalar x) {
-    /*
-     *
-     *     Psi (digamma) function (modified for Eigen)
-     *
-     *
-     * SYNOPSIS:
-     *
-     * double x, y, psi();
-     *
-     * y = psi( x );
-     *
-     *
-     * DESCRIPTION:
-     *
-     *              d      -
-     *   psi(x)  =  -- ln | (x)
-     *              dx
-     *
-     * is the logarithmic derivative of the gamma function.
-     * For integer x,
-     *                   n-1
-     *                    -
-     * psi(n) = -EUL  +   >  1/k.
-     *                    -
-     *                   k=1
-     *
-     * If x is negative, it is transformed to a positive argument by the
-     * reflection formula  psi(1-x) = psi(x) + pi cot(pi x).
-     * For general positive x, the argument is made greater than 10
-     * using the recurrence  psi(x+1) = psi(x) + 1/x.
-     * Then the following asymptotic expansion is applied:
-     *
-     *                           inf.   B
-     *                            -      2k
-     * psi(x) = log(x) - 1/2x -   >   -------
-     *                            -        2k
-     *                           k=1   2k x
-     *
-     * where the B2k are Bernoulli numbers.
-     *
-     * ACCURACY (float):
-     *    Relative error (except absolute when |psi| < 1):
-     * arithmetic   domain     # trials      peak         rms
-     *    IEEE      0,30        30000       1.3e-15     1.4e-16
-     *    IEEE      -30,0       40000       1.5e-15     2.2e-16
-     *
-     * ACCURACY (double):
-     *    Absolute error,  relative when |psi| > 1 :
-     * arithmetic   domain     # trials      peak         rms
-     *    IEEE      -33,0        30000      8.2e-7      1.2e-7
-     *    IEEE      0,33        100000      7.3e-7      7.7e-8
-     *
-     * ERROR MESSAGES:
-     *     message         condition      value returned
-     * psi singularity    x integer <=0      INFINITY
-     */
-
-    Scalar p, q, nz, s, w, y;
-    bool negative = false;
-
-    const Scalar maxnum = NumTraits<Scalar>::infinity();
-    const Scalar m_pi = Scalar(EIGEN_PI);
-
-    const Scalar zero = Scalar(0);
-    const Scalar one = Scalar(1);
-    const Scalar half = Scalar(0.5);
-    nz = zero;
-
-    if (x <= zero) {
-      negative = true;
-      q = x;
-      p = numext::floor(q);
-      if (p == q) {
-        return maxnum;
-      }
-      /* Remove the zeros of tan(m_pi x)
-       * by subtracting the nearest integer from x
-       */
-      nz = q - p;
-      if (nz != half) {
-        if (nz > half) {
-          p += one;
-          nz = q - p;
-        }
-        nz = m_pi / numext::tan(m_pi * nz);
-      }
-      else {
-        nz = zero;
-      }
-      x = one - x;
-    }
-
-    /* use the recurrence psi(x+1) = psi(x) + 1/x. */
-    s = x;
-    w = zero;
-    while (s < Scalar(10)) {
-      w += one / s;
-      s += one;
-    }
-
-    y = digamma_impl_maybe_poly<Scalar>::run(s);
-
-    y = numext::log(s) - (half / s) - y - w;
-
-    return (negative) ? y - nz : y;
-  }
-};
-
-/****************************************************************************
- * Implementation of erf, requires C++11/C99                                *
- ****************************************************************************/
-
-template <typename Scalar>
-struct erf_impl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-template <typename Scalar>
-struct erf_retval {
-  typedef Scalar type;
-};
-
-#if EIGEN_HAS_C99_MATH
-template <>
-struct erf_impl<float> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); }
-};
-
-template <>
-struct erf_impl<double> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double run(double x) { return ::erf(x); }
-};
-#endif  // EIGEN_HAS_C99_MATH
-
-/***************************************************************************
-* Implementation of erfc, requires C++11/C99                               *
-****************************************************************************/
-
-template <typename Scalar>
-struct erfc_impl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-template <typename Scalar>
-struct erfc_retval {
-  typedef Scalar type;
-};
-
-#if EIGEN_HAS_C99_MATH
-template <>
-struct erfc_impl<float> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); }
-};
-
-template <>
-struct erfc_impl<double> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double run(const double x) { return ::erfc(x); }
-};
-#endif  // EIGEN_HAS_C99_MATH
-
-/**************************************************************************************************************
- * Implementation of igammac (complemented incomplete gamma integral), based on Cephes but requires C++11/C99 *
- **************************************************************************************************************/
-
-template <typename Scalar>
-struct igammac_retval {
-  typedef Scalar type;
-};
-
-// NOTE: cephes_helper is also used to implement zeta
-template <typename Scalar>
-struct cephes_helper {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar machep() { assert(false && "machep not supported for this type"); return 0.0; }
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar big() { assert(false && "big not supported for this type"); return 0.0; }
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar biginv() { assert(false && "biginv not supported for this type"); return 0.0; }
-};
-
-template <>
-struct cephes_helper<float> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float machep() {
-    return NumTraits<float>::epsilon() / 2;  // 1.0 - machep == 1.0
-  }
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float big() {
-    // use epsneg (1.0 - epsneg == 1.0)
-    return 1.0f / (NumTraits<float>::epsilon() / 2);
-  }
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float biginv() {
-    // epsneg
-    return machep();
-  }
-};
-
-template <>
-struct cephes_helper<double> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double machep() {
-    return NumTraits<double>::epsilon() / 2;  // 1.0 - machep == 1.0
-  }
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double big() {
-    return 1.0 / NumTraits<double>::epsilon();
-  }
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double biginv() {
-    // inverse of eps
-    return NumTraits<double>::epsilon();
-  }
-};
-
-#if !EIGEN_HAS_C99_MATH
-
-template <typename Scalar>
-struct igammac_impl {
-  EIGEN_DEVICE_FUNC
-  static Scalar run(Scalar a, Scalar x) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-#else
-
-template <typename Scalar> struct igamma_impl;  // predeclare igamma_impl
-
-template <typename Scalar>
-struct igammac_impl {
-  EIGEN_DEVICE_FUNC
-  static Scalar run(Scalar a, Scalar x) {
-    /*  igamc()
-     *
-     *	Incomplete gamma integral (modified for Eigen)
-     *
-     *
-     *
-     * SYNOPSIS:
-     *
-     * double a, x, y, igamc();
-     *
-     * y = igamc( a, x );
-     *
-     * DESCRIPTION:
-     *
-     * The function is defined by
-     *
-     *
-     *  igamc(a,x)   =   1 - igam(a,x)
-     *
-     *                            inf.
-     *                              -
-     *                     1       | |  -t  a-1
-     *               =   -----     |   e   t   dt.
-     *                    -      | |
-     *                   | (a)    -
-     *                             x
-     *
-     *
-     * In this implementation both arguments must be positive.
-     * The integral is evaluated by either a power series or
-     * continued fraction expansion, depending on the relative
-     * values of a and x.
-     *
-     * ACCURACY (float):
-     *
-     *                      Relative error:
-     * arithmetic   domain     # trials      peak         rms
-     *    IEEE      0,30        30000       7.8e-6      5.9e-7
-     *
-     *
-     * ACCURACY (double):
-     *
-     * Tested at random a, x.
-     *                a         x                      Relative error:
-     * arithmetic   domain   domain     # trials      peak         rms
-     *    IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
-     *    IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
-     *
-     */
-    /*
-      Cephes Math Library Release 2.2: June, 1992
-      Copyright 1985, 1987, 1992 by Stephen L. Moshier
-      Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-    */
-    const Scalar zero = 0;
-    const Scalar one = 1;
-    const Scalar nan = NumTraits<Scalar>::quiet_NaN();
-
-    if ((x < zero) || (a <= zero)) {
-      // domain error
-      return nan;
-    }
-
-    if ((x < one) || (x < a)) {
-      /* The checks above ensure that we meet the preconditions for
-       * igamma_impl::Impl(), so call it, rather than igamma_impl::Run().
-       * Calling Run() would also work, but in that case the compiler may not be
-       * able to prove that igammac_impl::Run and igamma_impl::Run are not
-       * mutually recursive.  This leads to worse code, particularly on
-       * platforms like nvptx, where recursion is allowed only begrudgingly.
-       */
-      return (one - igamma_impl<Scalar>::Impl(a, x));
-    }
-
-    return Impl(a, x);
-  }
-
- private:
-  /* igamma_impl calls igammac_impl::Impl. */
-  friend struct igamma_impl<Scalar>;
-
-  /* Actually computes igamc(a, x).
-   *
-   * Preconditions:
-   *   a > 0
-   *   x >= 1
-   *   x >= a
-   */
-  EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) {
-    const Scalar zero = 0;
-    const Scalar one = 1;
-    const Scalar two = 2;
-    const Scalar machep = cephes_helper<Scalar>::machep();
-    const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
-    const Scalar big = cephes_helper<Scalar>::big();
-    const Scalar biginv = cephes_helper<Scalar>::biginv();
-    const Scalar inf = NumTraits<Scalar>::infinity();
-
-    Scalar ans, ax, c, yc, r, t, y, z;
-    Scalar pk, pkm1, pkm2, qk, qkm1, qkm2;
-
-    if (x == inf) return zero;  // std::isinf crashes on CUDA
-
-    /* Compute  x**a * exp(-x) / gamma(a)  */
-    ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
-    if (ax < -maxlog) {  // underflow
-      return zero;
-    }
-    ax = numext::exp(ax);
-
-    // continued fraction
-    y = one - a;
-    z = x + y + one;
-    c = zero;
-    pkm2 = one;
-    qkm2 = x;
-    pkm1 = x + one;
-    qkm1 = z * x;
-    ans = pkm1 / qkm1;
-
-    while (true) {
-      c += one;
-      y += one;
-      z += two;
-      yc = y * c;
-      pk = pkm1 * z - pkm2 * yc;
-      qk = qkm1 * z - qkm2 * yc;
-      if (qk != zero) {
-        r = pk / qk;
-        t = numext::abs((ans - r) / r);
-        ans = r;
-      } else {
-        t = one;
-      }
-      pkm2 = pkm1;
-      pkm1 = pk;
-      qkm2 = qkm1;
-      qkm1 = qk;
-      if (numext::abs(pk) > big) {
-        pkm2 *= biginv;
-        pkm1 *= biginv;
-        qkm2 *= biginv;
-        qkm1 *= biginv;
-      }
-      if (t <= machep) {
-        break;
-      }
-    }
-
-    return (ans * ax);
-  }
-};
-
-#endif  // EIGEN_HAS_C99_MATH
-
-/************************************************************************************************
- * Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 *
- ************************************************************************************************/
-
-template <typename Scalar>
-struct igamma_retval {
-  typedef Scalar type;
-};
-
-#if !EIGEN_HAS_C99_MATH
-
-template <typename Scalar>
-struct igamma_impl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-#else
-
-template <typename Scalar>
-struct igamma_impl {
-  EIGEN_DEVICE_FUNC
-  static Scalar run(Scalar a, Scalar x) {
-    /*	igam()
-     *	Incomplete gamma integral
-     *
-     *
-     *
-     * SYNOPSIS:
-     *
-     * double a, x, y, igam();
-     *
-     * y = igam( a, x );
-     *
-     * DESCRIPTION:
-     *
-     * The function is defined by
-     *
-     *                           x
-     *                            -
-     *                   1       | |  -t  a-1
-     *  igam(a,x)  =   -----     |   e   t   dt.
-     *                  -      | |
-     *                 | (a)    -
-     *                           0
-     *
-     *
-     * In this implementation both arguments must be positive.
-     * The integral is evaluated by either a power series or
-     * continued fraction expansion, depending on the relative
-     * values of a and x.
-     *
-     * ACCURACY (double):
-     *
-     *                      Relative error:
-     * arithmetic   domain     # trials      peak         rms
-     *    IEEE      0,30       200000       3.6e-14     2.9e-15
-     *    IEEE      0,100      300000       9.9e-14     1.5e-14
-     *
-     *
-     * ACCURACY (float):
-     *
-     *                      Relative error:
-     * arithmetic   domain     # trials      peak         rms
-     *    IEEE      0,30        20000       7.8e-6      5.9e-7
-     *
-     */
-    /*
-      Cephes Math Library Release 2.2: June, 1992
-      Copyright 1985, 1987, 1992 by Stephen L. Moshier
-      Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-    */
-
-
-    /* left tail of incomplete gamma function:
-     *
-     *          inf.      k
-     *   a  -x   -       x
-     *  x  e     >   ----------
-     *           -     -
-     *          k=0   | (a+k+1)
-     *
-     */
-    const Scalar zero = 0;
-    const Scalar one = 1;
-    const Scalar nan = NumTraits<Scalar>::quiet_NaN();
-
-    if (x == zero) return zero;
-
-    if ((x < zero) || (a <= zero)) {  // domain error
-      return nan;
-    }
-
-    if ((x > one) && (x > a)) {
-      /* The checks above ensure that we meet the preconditions for
-       * igammac_impl::Impl(), so call it, rather than igammac_impl::Run().
-       * Calling Run() would also work, but in that case the compiler may not be
-       * able to prove that igammac_impl::Run and igamma_impl::Run are not
-       * mutually recursive.  This leads to worse code, particularly on
-       * platforms like nvptx, where recursion is allowed only begrudgingly.
-       */
-      return (one - igammac_impl<Scalar>::Impl(a, x));
-    }
-
-    return Impl(a, x);
-  }
-
- private:
-  /* igammac_impl calls igamma_impl::Impl. */
-  friend struct igammac_impl<Scalar>;
-
-  /* Actually computes igam(a, x).
-   *
-   * Preconditions:
-   *   x > 0
-   *   a > 0
-   *   !(x > 1 && x > a)
-   */
-  EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) {
-    const Scalar zero = 0;
-    const Scalar one = 1;
-    const Scalar machep = cephes_helper<Scalar>::machep();
-    const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
-
-    Scalar ans, ax, c, r;
-
-    /* Compute  x**a * exp(-x) / gamma(a)  */
-    ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
-    if (ax < -maxlog) {
-      // underflow
-      return zero;
-    }
-    ax = numext::exp(ax);
-
-    /* power series */
-    r = a;
-    c = one;
-    ans = one;
-
-    while (true) {
-      r += one;
-      c *= x/r;
-      ans += c;
-      if (c/ans <= machep) {
-        break;
-      }
-    }
-
-    return (ans * ax / a);
-  }
-};
-
-#endif  // EIGEN_HAS_C99_MATH
-
-/*****************************************************************************
- * Implementation of Riemann zeta function of two arguments, based on Cephes *
- *****************************************************************************/
-
-template <typename Scalar>
-struct zeta_retval {
-    typedef Scalar type;
-};
-
-template <typename Scalar>
-struct zeta_impl_series {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-template <>
-struct zeta_impl_series<float> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE bool run(float& a, float& b, float& s, const float x, const float machep) {
-    int i = 0;  
-    while(i < 9)
-    {
-        i += 1;
-        a += 1.0f;
-        b = numext::pow( a, -x );
-        s += b;
-        if( numext::abs(b/s) < machep )
-            return true;
-    }
-    
-    //Return whether we are done
-    return false;
-  }
-};
-
-template <>
-struct zeta_impl_series<double> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE bool run(double& a, double& b, double& s, const double x, const double machep) {
-    int i = 0;  
-    while( (i < 9) || (a <= 9.0) )
-    {
-        i += 1;
-        a += 1.0;
-        b = numext::pow( a, -x );
-        s += b;
-        if( numext::abs(b/s) < machep )
-            return true;
-    }
-    
-    //Return whether we are done
-    return false;
-  }
-};
-    
-template <typename Scalar>
-struct zeta_impl {
-    EIGEN_DEVICE_FUNC
-    static Scalar run(Scalar x, Scalar q) {
-        /*							zeta.c
-         *
-         *	Riemann zeta function of two arguments
-         *
-         *
-         *
-         * SYNOPSIS:
-         *
-         * double x, q, y, zeta();
-         *
-         * y = zeta( x, q );
-         *
-         *
-         *
-         * DESCRIPTION:
-         *
-         *
-         *
-         *                 inf.
-         *                  -        -x
-         *   zeta(x,q)  =   >   (k+q)
-         *                  -
-         *                 k=0
-         *
-         * where x > 1 and q is not a negative integer or zero.
-         * The Euler-Maclaurin summation formula is used to obtain
-         * the expansion
-         *
-         *                n
-         *                -       -x
-         * zeta(x,q)  =   >  (k+q)
-         *                -
-         *               k=1
-         *
-         *           1-x                 inf.  B   x(x+1)...(x+2j)
-         *      (n+q)           1         -     2j
-         *  +  ---------  -  -------  +   >    --------------------
-         *        x-1              x      -                   x+2j+1
-         *                   2(n+q)      j=1       (2j)! (n+q)
-         *
-         * where the B2j are Bernoulli numbers.  Note that (see zetac.c)
-         * zeta(x,1) = zetac(x) + 1.
-         *
-         *
-         *
-         * ACCURACY:
-         *
-         * Relative error for single precision:
-         * arithmetic   domain     # trials      peak         rms
-         *    IEEE      0,25        10000       6.9e-7      1.0e-7
-         *
-         * Large arguments may produce underflow in powf(), in which
-         * case the results are inaccurate.
-         *
-         * REFERENCE:
-         *
-         * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
-         * Series, and Products, p. 1073; Academic Press, 1980.
-         *
-         */
-        
-        int i;
-        Scalar p, r, a, b, k, s, t, w;
-        
-        const Scalar A[] = {
-            Scalar(12.0),
-            Scalar(-720.0),
-            Scalar(30240.0),
-            Scalar(-1209600.0),
-            Scalar(47900160.0),
-            Scalar(-1.8924375803183791606e9), /*1.307674368e12/691*/
-            Scalar(7.47242496e10),
-            Scalar(-2.950130727918164224e12), /*1.067062284288e16/3617*/
-            Scalar(1.1646782814350067249e14), /*5.109094217170944e18/43867*/
-            Scalar(-4.5979787224074726105e15), /*8.028576626982912e20/174611*/
-            Scalar(1.8152105401943546773e17), /*1.5511210043330985984e23/854513*/
-            Scalar(-7.1661652561756670113e18) /*1.6938241367317436694528e27/236364091*/
-            };
-            
-        const Scalar maxnum = NumTraits<Scalar>::infinity();
-        const Scalar zero = 0.0, half = 0.5, one = 1.0;
-        const Scalar machep = cephes_helper<Scalar>::machep();
-        const Scalar nan = NumTraits<Scalar>::quiet_NaN();
-        
-        if( x == one )
-            return maxnum;
-        
-        if( x < one )
-        {
-            return nan;
-        }
-        
-        if( q <= zero )
-        {
-            if(q == numext::floor(q))
-            {
-                return maxnum;
-            }
-            p = x;
-            r = numext::floor(p);
-            if (p != r)
-                return nan;
-        }
-        
-        /* Permit negative q but continue sum until n+q > +9 .
-         * This case should be handled by a reflection formula.
-         * If q<0 and x is an integer, there is a relation to
-         * the polygamma function.
-         */
-        s = numext::pow( q, -x );
-        a = q;
-        b = zero;
-        // Run the summation in a helper function that is specific to the floating precision
-        if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) {
-            return s;
-        }
-        
-        w = a;
-        s += b*w/(x-one);
-        s -= half * b;
-        a = one;
-        k = zero;
-        for( i=0; i<12; i++ )
-        {
-            a *= x + k;
-            b /= w;
-            t = a*b/A[i];
-            s = s + t;
-            t = numext::abs(t/s);
-            if( t < machep ) {
-              break;
-            }
-            k += one;
-            a *= x + k;
-            b /= w;
-            k += one;
-        }
-        return s;
-  }
-};
-
-/****************************************************************************
- * Implementation of polygamma function, requires C++11/C99                 *
- ****************************************************************************/
-
-template <typename Scalar>
-struct polygamma_retval {
-    typedef Scalar type;
-};
-    
-#if !EIGEN_HAS_C99_MATH
-    
-template <typename Scalar>
-struct polygamma_impl {
-    EIGEN_DEVICE_FUNC
-    static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) {
-        EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                            THIS_TYPE_IS_NOT_SUPPORTED);
-        return Scalar(0);
-    }
-};
-    
-#else
-    
-template <typename Scalar>
-struct polygamma_impl {
-    EIGEN_DEVICE_FUNC
-    static Scalar run(Scalar n, Scalar x) {
-        Scalar zero = 0.0, one = 1.0;
-        Scalar nplus = n + one;
-        const Scalar nan = NumTraits<Scalar>::quiet_NaN();
-        
-        // Check that n is an integer
-        if (numext::floor(n) != n) {
-            return nan;
-        }
-        // Just return the digamma function for n = 1
-        else if (n == zero) {
-            return digamma_impl<Scalar>::run(x);
-        }
-        // Use the same implementation as scipy
-        else {
-            Scalar factorial = numext::exp(lgamma_impl<Scalar>::run(nplus));
-            return numext::pow(-one, nplus) * factorial * zeta_impl<Scalar>::run(nplus, x);
-        }
-  }
-};
-    
-#endif  // EIGEN_HAS_C99_MATH
-
-/************************************************************************************************
- * Implementation of betainc (incomplete beta integral), based on Cephes but requires C++11/C99 *
- ************************************************************************************************/
-
-template <typename Scalar>
-struct betainc_retval {
-  typedef Scalar type;
-};
-
-#if !EIGEN_HAS_C99_MATH
-
-template <typename Scalar>
-struct betainc_impl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-#else
-
-template <typename Scalar>
-struct betainc_impl {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(Scalar, Scalar, Scalar) {
-    /*	betaincf.c
-     *
-     *	Incomplete beta integral
-     *
-     *
-     * SYNOPSIS:
-     *
-     * float a, b, x, y, betaincf();
-     *
-     * y = betaincf( a, b, x );
-     *
-     *
-     * DESCRIPTION:
-     *
-     * Returns incomplete beta integral of the arguments, evaluated
-     * from zero to x.  The function is defined as
-     *
-     *                  x
-     *     -            -
-     *    | (a+b)      | |  a-1     b-1
-     *  -----------    |   t   (1-t)   dt.
-     *   -     -     | |
-     *  | (a) | (b)   -
-     *                 0
-     *
-     * The domain of definition is 0 <= x <= 1.  In this
-     * implementation a and b are restricted to positive values.
-     * The integral from x to 1 may be obtained by the symmetry
-     * relation
-     *
-     *    1 - betainc( a, b, x )  =  betainc( b, a, 1-x ).
-     *
-     * The integral is evaluated by a continued fraction expansion.
-     * If a < 1, the function calls itself recursively after a
-     * transformation to increase a to a+1.
-     *
-     * ACCURACY (float):
-     *
-     * Tested at random points (a,b,x) with a and b in the indicated
-     * interval and x between 0 and 1.
-     *
-     * arithmetic   domain     # trials      peak         rms
-     * Relative error:
-     *    IEEE       0,30       10000       3.7e-5      5.1e-6
-     *    IEEE       0,100      10000       1.7e-4      2.5e-5
-     * The useful domain for relative error is limited by underflow
-     * of the single precision exponential function.
-     * Absolute error:
-     *    IEEE       0,30      100000       2.2e-5      9.6e-7
-     *    IEEE       0,100      10000       6.5e-5      3.7e-6
-     *
-     * Larger errors may occur for extreme ratios of a and b.
-     *
-     * ACCURACY (double):
-     * arithmetic   domain     # trials      peak         rms
-     *    IEEE      0,5         10000       6.9e-15     4.5e-16
-     *    IEEE      0,85       250000       2.2e-13     1.7e-14
-     *    IEEE      0,1000      30000       5.3e-12     6.3e-13
-     *    IEEE      0,10000    250000       9.3e-11     7.1e-12
-     *    IEEE      0,100000    10000       8.7e-10     4.8e-11
-     * Outputs smaller than the IEEE gradual underflow threshold
-     * were excluded from these statistics.
-     *
-     * ERROR MESSAGES:
-     *   message         condition      value returned
-     * incbet domain      x<0, x>1          nan
-     * incbet underflow                     nan
-     */
-
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    return Scalar(0);
-  }
-};
-
-/* Continued fraction expansion #1 for incomplete beta integral (small_branch = True)
- * Continued fraction expansion #2 for incomplete beta integral (small_branch = False)
- */
-template <typename Scalar>
-struct incbeta_cfe {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x, bool small_branch) {
-    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, float>::value ||
-                         internal::is_same<Scalar, double>::value),
-                        THIS_TYPE_IS_NOT_SUPPORTED);
-    const Scalar big = cephes_helper<Scalar>::big();
-    const Scalar machep = cephes_helper<Scalar>::machep();
-    const Scalar biginv = cephes_helper<Scalar>::biginv();
-
-    const Scalar zero = 0;
-    const Scalar one = 1;
-    const Scalar two = 2;
-
-    Scalar xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-    Scalar k1, k2, k3, k4, k5, k6, k7, k8, k26update;
-    Scalar ans;
-    int n;
-
-    const int num_iters = (internal::is_same<Scalar, float>::value) ? 100 : 300;
-    const Scalar thresh =
-        (internal::is_same<Scalar, float>::value) ? machep : Scalar(3) * machep;
-    Scalar r = (internal::is_same<Scalar, float>::value) ? zero : one;
-
-    if (small_branch) {
-      k1 = a;
-      k2 = a + b;
-      k3 = a;
-      k4 = a + one;
-      k5 = one;
-      k6 = b - one;
-      k7 = k4;
-      k8 = a + two;
-      k26update = one;
-    } else {
-      k1 = a;
-      k2 = b - one;
-      k3 = a;
-      k4 = a + one;
-      k5 = one;
-      k6 = a + b;
-      k7 = a + one;
-      k8 = a + two;
-      k26update = -one;
-      x = x / (one - x);
-    }
-
-    pkm2 = zero;
-    qkm2 = one;
-    pkm1 = one;
-    qkm1 = one;
-    ans = one;
-    n = 0;
-
-    do {
-      xk = -(x * k1 * k2) / (k3 * k4);
-      pk = pkm1 + pkm2 * xk;
-      qk = qkm1 + qkm2 * xk;
-      pkm2 = pkm1;
-      pkm1 = pk;
-      qkm2 = qkm1;
-      qkm1 = qk;
-
-      xk = (x * k5 * k6) / (k7 * k8);
-      pk = pkm1 + pkm2 * xk;
-      qk = qkm1 + qkm2 * xk;
-      pkm2 = pkm1;
-      pkm1 = pk;
-      qkm2 = qkm1;
-      qkm1 = qk;
-
-      if (qk != zero) {
-        r = pk / qk;
-        if (numext::abs(ans - r) < numext::abs(r) * thresh) {
-          return r;
-        }
-        ans = r;
-      }
-
-      k1 += one;
-      k2 += k26update;
-      k3 += two;
-      k4 += two;
-      k5 += one;
-      k6 -= k26update;
-      k7 += two;
-      k8 += two;
-
-      if ((numext::abs(qk) + numext::abs(pk)) > big) {
-        pkm2 *= biginv;
-        pkm1 *= biginv;
-        qkm2 *= biginv;
-        qkm1 *= biginv;
-      }
-      if ((numext::abs(qk) < biginv) || (numext::abs(pk) < biginv)) {
-        pkm2 *= big;
-        pkm1 *= big;
-        qkm2 *= big;
-        qkm1 *= big;
-      }
-    } while (++n < num_iters);
-
-    return ans;
-  }
-};
-
-/* Helper functions depending on the Scalar type */
-template <typename Scalar>
-struct betainc_helper {};
-
-template <>
-struct betainc_helper<float> {
-  /* Core implementation, assumes a large (> 1.0) */
-  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE float incbsa(float aa, float bb,
-                                                            float xx) {
-    float ans, a, b, t, x, onemx;
-    bool reversed_a_b = false;
-
-    onemx = 1.0f - xx;
-
-    /* see if x is greater than the mean */
-    if (xx > (aa / (aa + bb))) {
-      reversed_a_b = true;
-      a = bb;
-      b = aa;
-      t = xx;
-      x = onemx;
-    } else {
-      a = aa;
-      b = bb;
-      t = onemx;
-      x = xx;
-    }
-
-    /* Choose expansion for optimal convergence */
-    if (b > 10.0f) {
-      if (numext::abs(b * x / a) < 0.3f) {
-        t = betainc_helper<float>::incbps(a, b, x);
-        if (reversed_a_b) t = 1.0f - t;
-        return t;
-      }
-    }
-
-    ans = x * (a + b - 2.0f) / (a - 1.0f);
-    if (ans < 1.0f) {
-      ans = incbeta_cfe<float>::run(a, b, x, true /* small_branch */);
-      t = b * numext::log(t);
-    } else {
-      ans = incbeta_cfe<float>::run(a, b, x, false /* small_branch */);
-      t = (b - 1.0f) * numext::log(t);
-    }
-
-    t += a * numext::log(x) + lgamma_impl<float>::run(a + b) -
-         lgamma_impl<float>::run(a) - lgamma_impl<float>::run(b);
-    t += numext::log(ans / a);
-    t = numext::exp(t);
-
-    if (reversed_a_b) t = 1.0f - t;
-    return t;
-  }
-
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE float incbps(float a, float b, float x) {
-    float t, u, y, s;
-    const float machep = cephes_helper<float>::machep();
-
-    y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a);
-    y -= lgamma_impl<float>::run(a) + lgamma_impl<float>::run(b);
-    y += lgamma_impl<float>::run(a + b);
-
-    t = x / (1.0f - x);
-    s = 0.0f;
-    u = 1.0f;
-    do {
-      b -= 1.0f;
-      if (b == 0.0f) {
-        break;
-      }
-      a += 1.0f;
-      u *= t * b / a;
-      s += u;
-    } while (numext::abs(u) > machep);
-
-    return numext::exp(y) * (1.0f + s);
-  }
-};
-
-template <>
-struct betainc_impl<float> {
-  EIGEN_DEVICE_FUNC
-  static float run(float a, float b, float x) {
-    const float nan = NumTraits<float>::quiet_NaN();
-    float ans, t;
-
-    if (a <= 0.0f) return nan;
-    if (b <= 0.0f) return nan;
-    if ((x <= 0.0f) || (x >= 1.0f)) {
-      if (x == 0.0f) return 0.0f;
-      if (x == 1.0f) return 1.0f;
-      // mtherr("betaincf", DOMAIN);
-      return nan;
-    }
-
-    /* transformation for small aa */
-    if (a <= 1.0f) {
-      ans = betainc_helper<float>::incbsa(a + 1.0f, b, x);
-      t = a * numext::log(x) + b * numext::log1p(-x) +
-          lgamma_impl<float>::run(a + b) - lgamma_impl<float>::run(a + 1.0f) -
-          lgamma_impl<float>::run(b);
-      return (ans + numext::exp(t));
-    } else {
-      return betainc_helper<float>::incbsa(a, b, x);
-    }
-  }
-};
-
-template <>
-struct betainc_helper<double> {
-  EIGEN_DEVICE_FUNC
-  static EIGEN_STRONG_INLINE double incbps(double a, double b, double x) {
-    const double machep = cephes_helper<double>::machep();
-
-    double s, t, u, v, n, t1, z, ai;
-
-    ai = 1.0 / a;
-    u = (1.0 - b) * x;
-    v = u / (a + 1.0);
-    t1 = v;
-    t = u;
-    n = 2.0;
-    s = 0.0;
-    z = machep * ai;
-    while (numext::abs(v) > z) {
-      u = (n - b) * x / n;
-      t *= u;
-      v = t / (a + n);
-      s += v;
-      n += 1.0;
-    }
-    s += t1;
-    s += ai;
-
-    u = a * numext::log(x);
-    // TODO: gamma() is not directly implemented in Eigen.
-    /*
-    if ((a + b) < maxgam && numext::abs(u) < maxlog) {
-      t = gamma(a + b) / (gamma(a) * gamma(b));
-      s = s * t * pow(x, a);
-    } else {
-    */
-    t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
-        lgamma_impl<double>::run(b) + u + numext::log(s);
-    return s = exp(t);
-  }
-};
-
-template <>
-struct betainc_impl<double> {
-  EIGEN_DEVICE_FUNC
-  static double run(double aa, double bb, double xx) {
-    const double nan = NumTraits<double>::quiet_NaN();
-    const double machep = cephes_helper<double>::machep();
-    // const double maxgam = 171.624376956302725;
-
-    double a, b, t, x, xc, w, y;
-    bool reversed_a_b = false;
-
-    if (aa <= 0.0 || bb <= 0.0) {
-      return nan;  // goto domerr;
-    }
-
-    if ((xx <= 0.0) || (xx >= 1.0)) {
-      if (xx == 0.0) return (0.0);
-      if (xx == 1.0) return (1.0);
-      // mtherr("incbet", DOMAIN);
-      return nan;
-    }
-
-    if ((bb * xx) <= 1.0 && xx <= 0.95) {
-      return betainc_helper<double>::incbps(aa, bb, xx);
-    }
-
-    w = 1.0 - xx;
-
-    /* Reverse a and b if x is greater than the mean. */
-    if (xx > (aa / (aa + bb))) {
-      reversed_a_b = true;
-      a = bb;
-      b = aa;
-      xc = xx;
-      x = w;
-    } else {
-      a = aa;
-      b = bb;
-      xc = w;
-      x = xx;
-    }
-
-    if (reversed_a_b && (b * x) <= 1.0 && x <= 0.95) {
-      t = betainc_helper<double>::incbps(a, b, x);
-      if (t <= machep) {
-        t = 1.0 - machep;
-      } else {
-        t = 1.0 - t;
-      }
-      return t;
-    }
-
-    /* Choose expansion for better convergence. */
-    y = x * (a + b - 2.0) - (a - 1.0);
-    if (y < 0.0) {
-      w = incbeta_cfe<double>::run(a, b, x, true /* small_branch */);
-    } else {
-      w = incbeta_cfe<double>::run(a, b, x, false /* small_branch */) / xc;
-    }
-
-    /* Multiply w by the factor
-         a      b   _             _     _
-        x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */
-
-    y = a * numext::log(x);
-    t = b * numext::log(xc);
-    // TODO: gamma is not directly implemented in Eigen.
-    /*
-    if ((a + b) < maxgam && numext::abs(y) < maxlog && numext::abs(t) < maxlog)
-    {
-      t = pow(xc, b);
-      t *= pow(x, a);
-      t /= a;
-      t *= w;
-      t *= gamma(a + b) / (gamma(a) * gamma(b));
-    } else {
-    */
-    /* Resort to logarithms.  */
-    y += t + lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
-         lgamma_impl<double>::run(b);
-    y += numext::log(w / a);
-    t = numext::exp(y);
-
-    /* } */
-    // done:
-
-    if (reversed_a_b) {
-      if (t <= machep) {
-        t = 1.0 - machep;
-      } else {
-        t = 1.0 - t;
-      }
-    }
-    return t;
-  }
-};
-
-#endif  // EIGEN_HAS_C99_MATH
-
-}  // end namespace internal
-
-namespace numext {
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar)
-    lgamma(const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar)
-    digamma(const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(zeta, Scalar)
-zeta(const Scalar& x, const Scalar& q) {
-    return EIGEN_MATHFUNC_IMPL(zeta, Scalar)::run(x, q);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(polygamma, Scalar)
-polygamma(const Scalar& n, const Scalar& x) {
-    return EIGEN_MATHFUNC_IMPL(polygamma, Scalar)::run(n, x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar)
-    erf(const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar)
-    erfc(const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar)
-    igamma(const Scalar& a, const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar)
-    igammac(const Scalar& a, const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(betainc, Scalar)
-    betainc(const Scalar& a, const Scalar& b, const Scalar& x) {
-  return EIGEN_MATHFUNC_IMPL(betainc, Scalar)::run(a, b, x);
-}
-
-}  // end namespace numext
-
-
-}  // end namespace Eigen
-
-#endif  // EIGEN_SPECIAL_FUNCTIONS_H
diff --git a/Eigen/src/Core/arch/CUDA/Half.h b/Eigen/src/Core/arch/CUDA/Half.h
index 87bdbfd1e..61692b83d 100644
--- a/Eigen/src/Core/arch/CUDA/Half.h
+++ b/Eigen/src/Core/arch/CUDA/Half.h
@@ -454,36 +454,6 @@ template <> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half maxi(const Eigen::
   return f1 < f2 ? b : a;
 #endif
 }
-
-#if EIGEN_HAS_C99_MATH
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half lgamma(const Eigen::half& a) {
-  return Eigen::half(Eigen::numext::lgamma(static_cast<float>(a)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half digamma(const Eigen::half& a) {
-  return Eigen::half(Eigen::numext::digamma(static_cast<float>(a)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half zeta(const Eigen::half& x, const Eigen::half& q) {
-  return Eigen::half(Eigen::numext::zeta(static_cast<float>(x), static_cast<float>(q)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half polygamma(const Eigen::half& n, const Eigen::half& x) {
-  return Eigen::half(Eigen::numext::polygamma(static_cast<float>(n), static_cast<float>(x)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erf(const Eigen::half& a) {
-  return Eigen::half(Eigen::numext::erf(static_cast<float>(a)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::half& a) {
-  return Eigen::half(Eigen::numext::erfc(static_cast<float>(a)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma(const Eigen::half& a, const Eigen::half& x) {
-  return Eigen::half(Eigen::numext::igamma(static_cast<float>(a), static_cast<float>(x)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igammac(const Eigen::half& a, const Eigen::half& x) {
-  return Eigen::half(Eigen::numext::igammac(static_cast<float>(a), static_cast<float>(x)));
-}
-template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half betainc(const Eigen::half& a, const Eigen::half& b, const Eigen::half& x) {
-  return Eigen::half(Eigen::numext::betainc(static_cast<float>(a), static_cast<float>(b), static_cast<float>(x)));
-}
-#endif
 } // end namespace numext
 
 } // end namespace Eigen
diff --git a/Eigen/src/Core/functors/BinaryFunctors.h b/Eigen/src/Core/functors/BinaryFunctors.h
index 2c1331208..dc3690444 100644
--- a/Eigen/src/Core/functors/BinaryFunctors.h
+++ b/Eigen/src/Core/functors/BinaryFunctors.h
@@ -429,57 +429,6 @@ template<> struct functor_traits<scalar_boolean_xor_op> {
   };
 };
 
-/** \internal
-  * \brief Template functor to compute the incomplete gamma function igamma(a, x)
-  *
-  * \sa class CwiseBinaryOp, Cwise::igamma
-  */
-template<typename Scalar> struct scalar_igamma_op : binary_op_base<Scalar,Scalar>
-{
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_op)
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const {
-    using numext::igamma; return igamma(a, x);
-  }
-  template<typename Packet>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const {
-    return internal::pigamma(a, x);
-  }
-};
-template<typename Scalar>
-struct functor_traits<scalar_igamma_op<Scalar> > {
-  enum {
-    // Guesstimate
-    Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasIGamma
-  };
-};
-
-
-/** \internal
-  * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x)
-  *
-  * \sa class CwiseBinaryOp, Cwise::igammac
-  */
-template<typename Scalar> struct scalar_igammac_op : binary_op_base<Scalar,Scalar>
-{
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_igammac_op)
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const {
-    using numext::igammac; return igammac(a, x);
-  }
-  template<typename Packet>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const
-  {
-    return internal::pigammac(a, x);
-  }
-};
-template<typename Scalar>
-struct functor_traits<scalar_igammac_op<Scalar> > {
-  enum {
-    // Guesstimate
-    Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasIGammac
-  };
-};
 
 
 //---------- binary functors bound to a constant, thus appearing as a unary functor ----------
diff --git a/Eigen/src/Core/functors/TernaryFunctors.h b/Eigen/src/Core/functors/TernaryFunctors.h
index 8b9e53062..b254e96c6 100644
--- a/Eigen/src/Core/functors/TernaryFunctors.h
+++ b/Eigen/src/Core/functors/TernaryFunctors.h
@@ -16,29 +16,7 @@ namespace internal {
 
 //---------- associative ternary functors ----------
 
-/** \internal
-  * \brief Template functor to compute the incomplete beta integral betainc(a, b, x)
-  *
-  */
-template<typename Scalar> struct scalar_betainc_op {
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_betainc_op)
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& x, const Scalar& a, const Scalar& b) const {
-    using numext::betainc; return betainc(x, a, b);
-  }
-  template<typename Packet>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& x, const Packet& a, const Packet& b) const
-  {
-    return internal::pbetainc(x, a, b);
-  }
-};
-template<typename Scalar>
-struct functor_traits<scalar_betainc_op<Scalar> > {
-  enum {
-    // Guesstimate
-    Cost = 400 * NumTraits<Scalar>::MulCost + 400 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasBetaInc
-  };
-};
+
 
 } // end namespace internal
 
diff --git a/Eigen/src/Core/functors/UnaryFunctors.h b/Eigen/src/Core/functors/UnaryFunctors.h
index a7d8c3b52..04208c9fe 100644
--- a/Eigen/src/Core/functors/UnaryFunctors.h
+++ b/Eigen/src/Core/functors/UnaryFunctors.h
@@ -473,142 +473,6 @@ struct functor_traits<scalar_asin_op<Scalar> >
 
 
 /** \internal
- * \brief Template functor to compute the natural log of the absolute
- * value of Gamma of a scalar
- * \sa class CwiseUnaryOp, Cwise::lgamma()
- */
-template<typename Scalar> struct scalar_lgamma_op {
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_lgamma_op)
-  EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const {
-    using numext::lgamma; return lgamma(a);
-  }
-  typedef typename packet_traits<Scalar>::type Packet;
-  EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::plgamma(a); }
-};
-template<typename Scalar>
-struct functor_traits<scalar_lgamma_op<Scalar> >
-{
-  enum {
-    // Guesstimate
-    Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasLGamma
-  };
-};
-
-/** \internal
- * \brief Template functor to compute psi, the derivative of lgamma of a scalar.
- * \sa class CwiseUnaryOp, Cwise::digamma()
- */
-template<typename Scalar> struct scalar_digamma_op {
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_digamma_op)
-  EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const {
-    using numext::digamma; return digamma(a);
-  }
-  typedef typename packet_traits<Scalar>::type Packet;
-  EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::pdigamma(a); }
-};
-template<typename Scalar>
-struct functor_traits<scalar_digamma_op<Scalar> >
-{
-  enum {
-    // Guesstimate
-    Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasDiGamma
-  };
-};
-    
-/** \internal
- * \brief Template functor to compute the Riemann Zeta function of two arguments.
- * \sa class CwiseUnaryOp, Cwise::zeta()
- */
-template<typename Scalar> struct scalar_zeta_op {
-    EIGEN_EMPTY_STRUCT_CTOR(scalar_zeta_op)
-    EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& x, const Scalar& q) const {
-        using numext::zeta; return zeta(x, q);
-    }
-    typedef typename packet_traits<Scalar>::type Packet;
-    EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); }
-};
-template<typename Scalar>
-struct functor_traits<scalar_zeta_op<Scalar> >
-{
-    enum {
-        // Guesstimate
-        Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
-        PacketAccess = packet_traits<Scalar>::HasZeta
-    };
-};
-
-/** \internal
- * \brief Template functor to compute the polygamma function.
- * \sa class CwiseUnaryOp, Cwise::polygamma()
- */
-template<typename Scalar> struct scalar_polygamma_op {
-    EIGEN_EMPTY_STRUCT_CTOR(scalar_polygamma_op)
-    EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& n, const Scalar& x) const {
-        using numext::polygamma; return polygamma(n, x);
-    }
-    typedef typename packet_traits<Scalar>::type Packet;
-    EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& n, const Packet& x) const { return internal::ppolygamma(n, x); }
-};
-template<typename Scalar>
-struct functor_traits<scalar_polygamma_op<Scalar> >
-{
-    enum {
-        // Guesstimate
-        Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
-        PacketAccess = packet_traits<Scalar>::HasPolygamma
-    };
-};
-
-/** \internal
- * \brief Template functor to compute the Gauss error function of a
- * scalar
- * \sa class CwiseUnaryOp, Cwise::erf()
- */
-template<typename Scalar> struct scalar_erf_op {
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_erf_op)
-  EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const {
-    using numext::erf; return erf(a);
-  }
-  typedef typename packet_traits<Scalar>::type Packet;
-  EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perf(a); }
-};
-template<typename Scalar>
-struct functor_traits<scalar_erf_op<Scalar> >
-{
-  enum {
-    // Guesstimate
-    Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasErf
-  };
-};
-
-/** \internal
- * \brief Template functor to compute the Complementary Error Function
- * of a scalar
- * \sa class CwiseUnaryOp, Cwise::erfc()
- */
-template<typename Scalar> struct scalar_erfc_op {
-  EIGEN_EMPTY_STRUCT_CTOR(scalar_erfc_op)
-  EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const {
-    using numext::erfc; return erfc(a);
-  }
-  typedef typename packet_traits<Scalar>::type Packet;
-  EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perfc(a); }
-};
-template<typename Scalar>
-struct functor_traits<scalar_erfc_op<Scalar> >
-{
-  enum {
-    // Guesstimate
-    Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
-    PacketAccess = packet_traits<Scalar>::HasErfc
-  };
-};
-
-
-/** \internal
   * \brief Template functor to compute the atan of a scalar
   * \sa class CwiseUnaryOp, ArrayBase::atan()
   */
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 1c90c0e2b..ea107393a 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -203,15 +203,21 @@ template<typename Scalar> struct scalar_random_op;
 template<typename Scalar> struct scalar_constant_op;
 template<typename Scalar> struct scalar_identity_op;
 template<typename Scalar,bool iscpx> struct scalar_sign_op;
-template<typename Scalar> struct scalar_igamma_op;
-template<typename Scalar> struct scalar_igammac_op;
-template<typename Scalar> struct scalar_betainc_op;
-
 template<typename Scalar,typename ScalarExponent> struct scalar_pow_op;
 template<typename LhsScalar,typename RhsScalar=LhsScalar> struct scalar_hypot_op;
 template<typename LhsScalar,typename RhsScalar=LhsScalar> struct scalar_product_op;
 template<typename LhsScalar,typename RhsScalar=LhsScalar> struct scalar_quotient_op;
 
+// SpecialFunctions module
+template<typename Scalar> struct scalar_lgamma_op;
+template<typename Scalar> struct scalar_digamma_op;
+template<typename Scalar> struct scalar_erf_op;
+template<typename Scalar> struct scalar_erfc_op;
+template<typename Scalar> struct scalar_igamma_op;
+template<typename Scalar> struct scalar_igammac_op;
+template<typename Scalar> struct scalar_zeta_op;
+template<typename Scalar> struct scalar_betainc_op;
+
 } // end namespace internal
 
 struct IOFormat;
diff --git a/Eigen/src/plugins/ArrayCwiseBinaryOps.h b/Eigen/src/plugins/ArrayCwiseBinaryOps.h
index 19e25ab62..62fb303d9 100644
--- a/Eigen/src/plugins/ArrayCwiseBinaryOps.h
+++ b/Eigen/src/plugins/ArrayCwiseBinaryOps.h
@@ -330,6 +330,8 @@ operator^(const EIGEN_CURRENT_STORAGE_BASE_CLASS<OtherDerived> &other) const
 #if 0
 /** \cpp11 \returns an expression of the coefficient-wise polygamma function.
   *
+  * \specialfunctions_module
+  *
   * It returns the \a n -th derivative of the digamma(psi) evaluated at \c *this.
   *
   * \warning Be careful with the order of the parameters: x.polygamma(n) is equivalent to polygamma(n,x)
@@ -346,6 +348,8 @@ polygamma(const EIGEN_CURRENT_STORAGE_BASE_CLASS<DerivedN> &n) const
 
 /** \returns an expression of the coefficient-wise zeta function.
   *
+  * \specialfunctions_module
+  *
   * It returns the Riemann zeta function of two arguments \c *this and \a q:
   *
   * \param *this is the exposent, it must be > 1
diff --git a/Eigen/src/plugins/ArrayCwiseUnaryOps.h b/Eigen/src/plugins/ArrayCwiseUnaryOps.h
index 9e42bb540..db02e299c 100644
--- a/Eigen/src/plugins/ArrayCwiseUnaryOps.h
+++ b/Eigen/src/plugins/ArrayCwiseUnaryOps.h
@@ -22,10 +22,6 @@ typedef CwiseUnaryOp<internal::scalar_atan_op<Scalar>, const Derived> AtanReturn
 typedef CwiseUnaryOp<internal::scalar_tanh_op<Scalar>, const Derived> TanhReturnType;
 typedef CwiseUnaryOp<internal::scalar_sinh_op<Scalar>, const Derived> SinhReturnType;
 typedef CwiseUnaryOp<internal::scalar_cosh_op<Scalar>, const Derived> CoshReturnType;
-typedef CwiseUnaryOp<internal::scalar_lgamma_op<Scalar>, const Derived> LgammaReturnType;
-typedef CwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived> DigammaReturnType;
-typedef CwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived> ErfReturnType;
-typedef CwiseUnaryOp<internal::scalar_erfc_op<Scalar>, const Derived> ErfcReturnType;
 typedef CwiseUnaryOp<internal::scalar_square_op<Scalar>, const Derived> SquareReturnType;
 typedef CwiseUnaryOp<internal::scalar_cube_op<Scalar>, const Derived> CubeReturnType;
 typedef CwiseUnaryOp<internal::scalar_round_op<Scalar>, const Derived> RoundReturnType;
@@ -324,77 +320,6 @@ cosh() const
   return CoshReturnType(derived());
 }
 
-/** \cpp11 \returns an expression of the coefficient-wise ln(|gamma(*this)|).
-  *
-  * Example: \include Cwise_lgamma.cpp
-  * Output: \verbinclude Cwise_lgamma.out
-  *
-  * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-  * or float/double in non c++11 mode, the user has to provide implementations of lgamma(T) for any scalar
-  * type T to be supported.
-  *
-  * \sa digamma()
-  */
-EIGEN_DEVICE_FUNC
-inline const LgammaReturnType
-lgamma() const
-{
-  return LgammaReturnType(derived());
-}
-
-/** \returns an expression of the coefficient-wise digamma (psi, derivative of lgamma).
-  *
-  * \note This function supports only float and double scalar types. To support other scalar types,
-  * the user has to provide implementations of digamma(T) for any scalar
-  * type T to be supported.
-  *
-  * \sa Eigen::digamma(), Eigen::polygamma(), lgamma()
-  */
-EIGEN_DEVICE_FUNC
-inline const DigammaReturnType
-digamma() const
-{
-  return DigammaReturnType(derived());
-}
-
-/** \cpp11 \returns an expression of the coefficient-wise Gauss error
-  * function of *this.
-  *
-  * Example: \include Cwise_erf.cpp
-  * Output: \verbinclude Cwise_erf.out
-  *
-  * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-  * or float/double in non c++11 mode, the user has to provide implementations of erf(T) for any scalar
-  * type T to be supported.
-  *
-  * \sa erfc()
-  */
-EIGEN_DEVICE_FUNC
-inline const ErfReturnType
-erf() const
-{
-  return ErfReturnType(derived());
-}
-
-/** \cpp11 \returns an expression of the coefficient-wise Complementary error
-  * function of *this.
-  *
-  * Example: \include Cwise_erfc.cpp
-  * Output: \verbinclude Cwise_erfc.out
-  *
-  * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
-  * or float/double in non c++11 mode, the user has to provide implementations of erfc(T) for any scalar
-  * type T to be supported.
-  *
-  * \sa erf()
-  */
-EIGEN_DEVICE_FUNC
-inline const ErfcReturnType
-erfc() const
-{
-  return ErfcReturnType(derived());
-}
-
 /** \returns an expression of the coefficient-wise inverse of *this.
   *
   * Example: \include Cwise_inverse.cpp
@@ -538,3 +463,90 @@ operator!() const
                       THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_OF_BOOL);
   return BooleanNotReturnType(derived());
 }
+
+
+// --- SpecialFunctions module ---
+
+typedef CwiseUnaryOp<internal::scalar_lgamma_op<Scalar>, const Derived> LgammaReturnType;
+typedef CwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived> DigammaReturnType;
+typedef CwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived> ErfReturnType;
+typedef CwiseUnaryOp<internal::scalar_erfc_op<Scalar>, const Derived> ErfcReturnType;
+
+/** \cpp11 \returns an expression of the coefficient-wise ln(|gamma(*this)|).
+  *
+  * \specialfunctions_module
+  *
+  * Example: \include Cwise_lgamma.cpp
+  * Output: \verbinclude Cwise_lgamma.out
+  *
+  * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
+  * or float/double in non c++11 mode, the user has to provide implementations of lgamma(T) for any scalar
+  * type T to be supported.
+  *
+  * \sa digamma()
+  */
+EIGEN_DEVICE_FUNC
+inline const LgammaReturnType
+lgamma() const
+{
+  return LgammaReturnType(derived());
+}
+
+/** \returns an expression of the coefficient-wise digamma (psi, derivative of lgamma).
+  *
+  * \specialfunctions_module
+  *
+  * \note This function supports only float and double scalar types. To support other scalar types,
+  * the user has to provide implementations of digamma(T) for any scalar
+  * type T to be supported.
+  *
+  * \sa Eigen::digamma(), Eigen::polygamma(), lgamma()
+  */
+EIGEN_DEVICE_FUNC
+inline const DigammaReturnType
+digamma() const
+{
+  return DigammaReturnType(derived());
+}
+
+/** \cpp11 \returns an expression of the coefficient-wise Gauss error
+  * function of *this.
+  *
+  * \specialfunctions_module
+  *
+  * Example: \include Cwise_erf.cpp
+  * Output: \verbinclude Cwise_erf.out
+  *
+  * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
+  * or float/double in non c++11 mode, the user has to provide implementations of erf(T) for any scalar
+  * type T to be supported.
+  *
+  * \sa erfc()
+  */
+EIGEN_DEVICE_FUNC
+inline const ErfReturnType
+erf() const
+{
+  return ErfReturnType(derived());
+}
+
+/** \cpp11 \returns an expression of the coefficient-wise Complementary error
+  * function of *this.
+  *
+  * \specialfunctions_module
+  *
+  * Example: \include Cwise_erfc.cpp
+  * Output: \verbinclude Cwise_erfc.out
+  *
+  * \note This function supports only float and double scalar types in c++11 mode. To support other scalar types,
+  * or float/double in non c++11 mode, the user has to provide implementations of erfc(T) for any scalar
+  * type T to be supported.
+  *
+  * \sa erf()
+  */
+EIGEN_DEVICE_FUNC
+inline const ErfcReturnType
+erfc() const
+{
+  return ErfcReturnType(derived());
+}