From 2f7e2614e773dde8a84156b4e3864474af8b53d6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 8 Jul 2016 11:13:55 +0200 Subject: bug #1232: refactor special functions as a new SpecialFunctions module, currently in unsupported/. --- Eigen/Core | 1 - Eigen/src/Core/GenericPacketMath.h | 36 - Eigen/src/Core/GlobalFunctions.h | 105 -- Eigen/src/Core/SpecialFunctions.h | 1551 -------------------- Eigen/src/Core/arch/CUDA/Half.h | 30 - Eigen/src/Core/functors/BinaryFunctors.h | 51 - Eigen/src/Core/functors/TernaryFunctors.h | 24 +- Eigen/src/Core/functors/UnaryFunctors.h | 136 -- Eigen/src/Core/util/ForwardDeclarations.h | 14 +- Eigen/src/plugins/ArrayCwiseBinaryOps.h | 4 + Eigen/src/plugins/ArrayCwiseUnaryOps.h | 162 +- doc/Doxyfile.in | 1 + test/array.cpp | 326 +--- unsupported/Eigen/CMakeLists.txt | 1 + unsupported/Eigen/CXX11/Tensor | 1 + unsupported/Eigen/SpecialFunctions | 57 + unsupported/Eigen/src/CMakeLists.txt | 1 + .../Eigen/src/SpecialFunctions/CMakeLists.txt | 6 + .../SpecialFunctions/SpecialFunctionsArrayAPI.h | 124 ++ .../SpecialFunctions/SpecialFunctionsFunctors.h | 236 +++ .../src/SpecialFunctions/SpecialFunctionsHalf.h | 47 + .../src/SpecialFunctions/SpecialFunctionsImpl.h | 1551 ++++++++++++++++++++ .../SpecialFunctions/SpecialFunctionsPacketMath.h | 58 + unsupported/test/CMakeLists.txt | 1 + unsupported/test/special_functions.cpp | 345 +++++ 25 files changed, 2532 insertions(+), 2337 deletions(-) delete mode 100644 Eigen/src/Core/SpecialFunctions.h create mode 100644 unsupported/Eigen/SpecialFunctions create mode 100644 unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt create mode 100644 unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h create mode 100644 unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h create mode 100644 unsupported/test/special_functions.cpp diff --git a/Eigen/Core b/Eigen/Core index 1f3a6504e..946ed0677 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -328,7 +328,6 @@ using std::ptrdiff_t; #include "src/Core/NumTraits.h" #include "src/Core/MathFunctions.h" -#include "src/Core/SpecialFunctions.h" #include "src/Core/GenericPacketMath.h" #if defined EIGEN_VECTORIZE_AVX 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 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 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 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 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 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 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 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 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 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 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 - inline const Eigen::CwiseBinaryOp, const Derived, const ExponentDerived> - igamma(const Eigen::ArrayBase& a, const Eigen::ArrayBase& x) - { - return Eigen::CwiseBinaryOp, 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 - inline const Eigen::CwiseBinaryOp, const Derived, const ExponentDerived> - igammac(const Eigen::ArrayBase& a, const Eigen::ArrayBase& x) - { - return Eigen::CwiseBinaryOp, 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 - inline const Eigen::CwiseBinaryOp, const DerivedN, const DerivedX> - polygamma(const Eigen::ArrayBase& n, const Eigen::ArrayBase& x) - { - return Eigen::CwiseBinaryOp, 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 - inline const Eigen::CwiseTernaryOp, const ArgADerived, const ArgBDerived, const ArgXDerived> - betainc(const Eigen::ArrayBase& a, const Eigen::ArrayBase& b, const Eigen::ArrayBase& x) - { - return Eigen::CwiseTernaryOp, 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 - inline const Eigen::CwiseBinaryOp, const DerivedX, const DerivedQ> - zeta(const Eigen::ArrayBase& x, const Eigen::ArrayBase& q) - { - return Eigen::CwiseBinaryOp, 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 -// -// 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( 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 -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::run(x, coef) * x + coef[N]; - } -}; - -template -struct polevl { - 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 -struct lgamma_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -template -struct lgamma_retval { - typedef Scalar type; -}; - -#if EIGEN_HAS_C99_MATH -template <> -struct lgamma_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float run(float x) { return ::lgammaf(x); } -}; - -template <> -struct lgamma_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double run(double x) { return ::lgamma(x); } -}; -#endif - -/**************************************************************************** - * Implementation of digamma (psi), based on Cephes * - ****************************************************************************/ - -template -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 -struct digamma_impl_maybe_poly { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - - -template <> -struct digamma_impl_maybe_poly { - 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::run(z, A); - } else return 0.0f; - } -}; - -template <> -struct digamma_impl_maybe_poly { - 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::run(z, A); - } - else return 0.0; - } -}; - -template -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::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::run(s); - - y = numext::log(s) - (half / s) - y - w; - - return (negative) ? y - nz : y; - } -}; - -/**************************************************************************** - * Implementation of erf, requires C++11/C99 * - ****************************************************************************/ - -template -struct erf_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -template -struct erf_retval { - typedef Scalar type; -}; - -#if EIGEN_HAS_C99_MATH -template <> -struct erf_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); } -}; - -template <> -struct erf_impl { - 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 -struct erfc_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -template -struct erfc_retval { - typedef Scalar type; -}; - -#if EIGEN_HAS_C99_MATH -template <> -struct erfc_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); } -}; - -template <> -struct erfc_impl { - 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 -struct igammac_retval { - typedef Scalar type; -}; - -// NOTE: cephes_helper is also used to implement zeta -template -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 { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float machep() { - return NumTraits::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::epsilon() / 2); - } - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float biginv() { - // epsneg - return machep(); - } -}; - -template <> -struct cephes_helper { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double machep() { - return NumTraits::epsilon() / 2; // 1.0 - machep == 1.0 - } - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double big() { - return 1.0 / NumTraits::epsilon(); - } - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double biginv() { - // inverse of eps - return NumTraits::epsilon(); - } -}; - -#if !EIGEN_HAS_C99_MATH - -template -struct igammac_impl { - EIGEN_DEVICE_FUNC - static Scalar run(Scalar a, Scalar x) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -#else - -template struct igamma_impl; // predeclare igamma_impl - -template -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::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::Impl(a, x)); - } - - return Impl(a, x); - } - - private: - /* igamma_impl calls igammac_impl::Impl. */ - friend struct igamma_impl; - - /* 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::machep(); - const Scalar maxlog = numext::log(NumTraits::highest()); - const Scalar big = cephes_helper::big(); - const Scalar biginv = cephes_helper::biginv(); - const Scalar inf = NumTraits::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::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 -struct igamma_retval { - typedef Scalar type; -}; - -#if !EIGEN_HAS_C99_MATH - -template -struct igamma_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -#else - -template -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::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::Impl(a, x)); - } - - return Impl(a, x); - } - - private: - /* igammac_impl calls igamma_impl::Impl. */ - friend struct igammac_impl; - - /* 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::machep(); - const Scalar maxlog = numext::log(NumTraits::highest()); - - Scalar ans, ax, c, r; - - /* Compute x**a * exp(-x) / gamma(a) */ - ax = a * numext::log(x) - x - lgamma_impl::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 -struct zeta_retval { - typedef Scalar type; -}; - -template -struct zeta_impl_series { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -template <> -struct zeta_impl_series { - 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 { - 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 -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::infinity(); - const Scalar zero = 0.0, half = 0.5, one = 1.0; - const Scalar machep = cephes_helper::machep(); - const Scalar nan = NumTraits::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::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 -struct polygamma_retval { - typedef Scalar type; -}; - -#if !EIGEN_HAS_C99_MATH - -template -struct polygamma_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -#else - -template -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::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::run(x); - } - // Use the same implementation as scipy - else { - Scalar factorial = numext::exp(lgamma_impl::run(nplus)); - return numext::pow(-one, nplus) * factorial * zeta_impl::run(nplus, x); - } - } -}; - -#endif // EIGEN_HAS_C99_MATH - -/************************************************************************************************ - * Implementation of betainc (incomplete beta integral), based on Cephes but requires C++11/C99 * - ************************************************************************************************/ - -template -struct betainc_retval { - typedef Scalar type; -}; - -#if !EIGEN_HAS_C99_MATH - -template -struct betainc_impl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) { - EIGEN_STATIC_ASSERT((internal::is_same::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); - } -}; - -#else - -template -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::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 -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::value || - internal::is_same::value), - THIS_TYPE_IS_NOT_SUPPORTED); - const Scalar big = cephes_helper::big(); - const Scalar machep = cephes_helper::machep(); - const Scalar biginv = cephes_helper::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::value) ? 100 : 300; - const Scalar thresh = - (internal::is_same::value) ? machep : Scalar(3) * machep; - Scalar r = (internal::is_same::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 -struct betainc_helper {}; - -template <> -struct betainc_helper { - /* 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::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::run(a, b, x, true /* small_branch */); - t = b * numext::log(t); - } else { - ans = incbeta_cfe::run(a, b, x, false /* small_branch */); - t = (b - 1.0f) * numext::log(t); - } - - t += a * numext::log(x) + lgamma_impl::run(a + b) - - lgamma_impl::run(a) - lgamma_impl::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::machep(); - - y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a); - y -= lgamma_impl::run(a) + lgamma_impl::run(b); - y += lgamma_impl::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 { - EIGEN_DEVICE_FUNC - static float run(float a, float b, float x) { - const float nan = NumTraits::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::incbsa(a + 1.0f, b, x); - t = a * numext::log(x) + b * numext::log1p(-x) + - lgamma_impl::run(a + b) - lgamma_impl::run(a + 1.0f) - - lgamma_impl::run(b); - return (ans + numext::exp(t)); - } else { - return betainc_helper::incbsa(a, b, x); - } - } -}; - -template <> -struct betainc_helper { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double incbps(double a, double b, double x) { - const double machep = cephes_helper::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::run(a + b) - lgamma_impl::run(a) - - lgamma_impl::run(b) + u + numext::log(s); - return s = exp(t); - } -}; - -template <> -struct betainc_impl { - EIGEN_DEVICE_FUNC - static double run(double aa, double bb, double xx) { - const double nan = NumTraits::quiet_NaN(); - const double machep = cephes_helper::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::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::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::run(a, b, x, true /* small_branch */); - } else { - w = incbeta_cfe::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::run(a + b) - lgamma_impl::run(a) - - lgamma_impl::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 -EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar) - lgamma(const Scalar& x) { - return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x); -} - -template -EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar) - digamma(const Scalar& x) { - return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x); -} - -template -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 -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 -EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar) - erf(const Scalar& x) { - return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x); -} - -template -EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) - erfc(const Scalar& x) { - return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x); -} - -template -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 -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 -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(a))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half digamma(const Eigen::half& a) { - return Eigen::half(Eigen::numext::digamma(static_cast(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(x), static_cast(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(n), static_cast(x))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erf(const Eigen::half& a) { - return Eigen::half(Eigen::numext::erf(static_cast(a))); -} -template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::half& a) { - return Eigen::half(Eigen::numext::erfc(static_cast(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(a), static_cast(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(a), static_cast(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(a), static_cast(b), static_cast(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 { }; }; -/** \internal - * \brief Template functor to compute the incomplete gamma function igamma(a, x) - * - * \sa class CwiseBinaryOp, Cwise::igamma - */ -template struct scalar_igamma_op : binary_op_base -{ - 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 - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const { - return internal::pigamma(a, x); - } -}; -template -struct functor_traits > { - enum { - // Guesstimate - Cost = 20 * NumTraits::MulCost + 10 * NumTraits::AddCost, - PacketAccess = packet_traits::HasIGamma - }; -}; - - -/** \internal - * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x) - * - * \sa class CwiseBinaryOp, Cwise::igammac - */ -template struct scalar_igammac_op : binary_op_base -{ - 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 - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const - { - return internal::pigammac(a, x); - } -}; -template -struct functor_traits > { - enum { - // Guesstimate - Cost = 20 * NumTraits::MulCost + 10 * NumTraits::AddCost, - PacketAccess = packet_traits::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 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 - 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 -struct functor_traits > { - enum { - // Guesstimate - Cost = 400 * NumTraits::MulCost + 400 * NumTraits::AddCost, - PacketAccess = packet_traits::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 @@ -472,142 +472,6 @@ struct functor_traits > }; -/** \internal - * \brief Template functor to compute the natural log of the absolute - * value of Gamma of a scalar - * \sa class CwiseUnaryOp, Cwise::lgamma() - */ -template 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::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::plgamma(a); } -}; -template -struct functor_traits > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, - PacketAccess = packet_traits::HasLGamma - }; -}; - -/** \internal - * \brief Template functor to compute psi, the derivative of lgamma of a scalar. - * \sa class CwiseUnaryOp, Cwise::digamma() - */ -template 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::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::pdigamma(a); } -}; -template -struct functor_traits > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, - PacketAccess = packet_traits::HasDiGamma - }; -}; - -/** \internal - * \brief Template functor to compute the Riemann Zeta function of two arguments. - * \sa class CwiseUnaryOp, Cwise::zeta() - */ -template 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::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); } -}; -template -struct functor_traits > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, - PacketAccess = packet_traits::HasZeta - }; -}; - -/** \internal - * \brief Template functor to compute the polygamma function. - * \sa class CwiseUnaryOp, Cwise::polygamma() - */ -template 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::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& n, const Packet& x) const { return internal::ppolygamma(n, x); } -}; -template -struct functor_traits > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, - PacketAccess = packet_traits::HasPolygamma - }; -}; - -/** \internal - * \brief Template functor to compute the Gauss error function of a - * scalar - * \sa class CwiseUnaryOp, Cwise::erf() - */ -template 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::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perf(a); } -}; -template -struct functor_traits > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, - PacketAccess = packet_traits::HasErf - }; -}; - -/** \internal - * \brief Template functor to compute the Complementary Error Function - * of a scalar - * \sa class CwiseUnaryOp, Cwise::erfc() - */ -template 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::type Packet; - EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perfc(a); } -}; -template -struct functor_traits > -{ - enum { - // Guesstimate - Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, - PacketAccess = packet_traits::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 struct scalar_random_op; template struct scalar_constant_op; template struct scalar_identity_op; template struct scalar_sign_op; -template struct scalar_igamma_op; -template struct scalar_igammac_op; -template struct scalar_betainc_op; - template struct scalar_pow_op; template struct scalar_hypot_op; template struct scalar_product_op; template struct scalar_quotient_op; +// SpecialFunctions module +template struct scalar_lgamma_op; +template struct scalar_digamma_op; +template struct scalar_erf_op; +template struct scalar_erfc_op; +template struct scalar_igamma_op; +template struct scalar_igammac_op; +template struct scalar_zeta_op; +template 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 @@ -329,6 +329,8 @@ operator^(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const // NOTE disabled until we agree on argument order #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. * @@ -345,6 +347,8 @@ polygamma(const EIGEN_CURRENT_STORAGE_BASE_CLASS &n) const #endif /** \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: * 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, const Derived> AtanReturn typedef CwiseUnaryOp, const Derived> TanhReturnType; typedef CwiseUnaryOp, const Derived> SinhReturnType; typedef CwiseUnaryOp, const Derived> CoshReturnType; -typedef CwiseUnaryOp, const Derived> LgammaReturnType; -typedef CwiseUnaryOp, const Derived> DigammaReturnType; -typedef CwiseUnaryOp, const Derived> ErfReturnType; -typedef CwiseUnaryOp, const Derived> ErfcReturnType; typedef CwiseUnaryOp, const Derived> SquareReturnType; typedef CwiseUnaryOp, const Derived> CubeReturnType; typedef CwiseUnaryOp, 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, const Derived> LgammaReturnType; +typedef CwiseUnaryOp, const Derived> DigammaReturnType; +typedef CwiseUnaryOp, const Derived> ErfReturnType; +typedef CwiseUnaryOp, 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()); +} diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 0c3673f89..cfc15a0f9 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -216,6 +216,7 @@ ALIASES = "only_for_vectors=This is only for vectors (either row- "lu_module=This is defined in the %LU module. \code #include \endcode" \ "qr_module=This is defined in the %QR module. \code #include \endcode" \ "svd_module=This is defined in the %SVD module. \code #include \endcode" \ + "specialfunctions_module=This is defined in the SpecialFunctions module. \code #include \endcode" \ "label=\bug" \ "matrixworld=*" \ "arrayworld=*" \ diff --git a/test/array.cpp b/test/array.cpp index 6347c8067..d734e604a 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -234,12 +234,7 @@ template void array_real(const ArrayType& m) VERIFY_IS_APPROX(m1.sinh(), sinh(m1)); VERIFY_IS_APPROX(m1.cosh(), cosh(m1)); VERIFY_IS_APPROX(m1.tanh(), tanh(m1)); -#if EIGEN_HAS_C99_MATH - VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1)); - VERIFY_IS_APPROX(m1.digamma(), digamma(m1)); - VERIFY_IS_APPROX(m1.erf(), erf(m1)); - VERIFY_IS_APPROX(m1.erfc(), erfc(m1)); -#endif // EIGEN_HAS_C99_MATH + VERIFY_IS_APPROX(m1.arg(), arg(m1)); VERIFY_IS_APPROX(m1.round(), round(m1)); VERIFY_IS_APPROX(m1.floor(), floor(m1)); @@ -313,88 +308,6 @@ template void array_real(const ArrayType& m) m1 += ArrayType::Constant(rows,cols,Scalar(tiny)); VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse()); - - -#if EIGEN_HAS_C99_MATH - // check special functions (comparing against numpy implementation) - if (!NumTraits::IsComplex) - { - - { - // Test various propreties of igamma & igammac. These are normalized - // gamma integrals where - // igammac(a, x) = Gamma(a, x) / Gamma(a) - // igamma(a, x) = gamma(a, x) / Gamma(a) - // where Gamma and gamma are considered the standard unnormalized - // upper and lower incomplete gamma functions, respectively. - ArrayType a = m1.abs() + 2; - ArrayType x = m2.abs() + 2; - ArrayType zero = ArrayType::Zero(rows, cols); - ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0)); - ArrayType a_m1 = a - one; - ArrayType Gamma_a_x = Eigen::igammac(a, x) * a.lgamma().exp(); - ArrayType Gamma_a_m1_x = Eigen::igammac(a_m1, x) * a_m1.lgamma().exp(); - ArrayType gamma_a_x = Eigen::igamma(a, x) * a.lgamma().exp(); - ArrayType gamma_a_m1_x = Eigen::igamma(a_m1, x) * a_m1.lgamma().exp(); - - // Gamma(a, 0) == Gamma(a) - VERIFY_IS_APPROX(Eigen::igammac(a, zero), one); - - // Gamma(a, x) + gamma(a, x) == Gamma(a) - VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp()); - - // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) - VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).exp()); - - // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) - VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).exp()); - } - - // Check exact values of igamma and igammac against a third party calculation. - Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; - Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; - - // location i*6+j corresponds to a_s[i], x_s[j]. - Scalar nan = std::numeric_limits::quiet_NaN(); - Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan}, - {0.0, 0.6321205588285578, 0.7768698398515702, - 0.9816843611112658, 9.999500016666262e-05, 1.0}, - {0.0, 0.4275932955291202, 0.608374823728911, - 0.9539882943107686, 7.522076445089201e-07, 1.0}, - {0.0, 0.01898815687615381, 0.06564245437845008, - 0.5665298796332909, 4.166333347221828e-18, 1.0}, - {0.0, 0.9999780593618628, 0.9999899967080838, - 0.9999996219837988, 0.9991370418689945, 1.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}}; - Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan}, - {1.0, 0.36787944117144233, 0.22313016014842982, - 0.018315638888734182, 0.9999000049998333, 0.0}, - {1.0, 0.5724067044708798, 0.3916251762710878, - 0.04601170568923136, 0.9999992477923555, 0.0}, - {1.0, 0.9810118431238462, 0.9343575456215499, - 0.4334701203667089, 1.0, 0.0}, - {1.0, 2.1940638138146658e-05, 1.0003291916285e-05, - 3.7801620118431334e-07, 0.0008629581310054535, - 0.0}, - {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}}; - for (int i = 0; i < 6; ++i) { - for (int j = 0; j < 6; ++j) { - if ((std::isnan)(igamma_s[i][j])) { - VERIFY((std::isnan)(numext::igamma(a_s[i], x_s[j]))); - } else { - VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]); - } - - if ((std::isnan)(igammac_s[i][j])) { - VERIFY((std::isnan)(numext::igammac(a_s[i], x_s[j]))); - } else { - VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]); - } - } - } - } -#endif // EIGEN_HAS_C99_MATH - // check inplace transpose m3 = m1; m3.transposeInPlace(); @@ -537,242 +450,8 @@ template void min_max(const ArrayType& m) } -template -void verify_component_wise(const X& x, const Y& y) -{ - for(Index i=0; i void array_special_functions() -{ - using std::abs; - using std::sqrt; - typedef typename ArrayType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - - Scalar plusinf = std::numeric_limits::infinity(); - Scalar nan = std::numeric_limits::quiet_NaN(); - - // Check the zeta function against scipy.special.zeta - { - ArrayType x(7), q(7), res(7), ref(7); - x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9; - q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345; - ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan; - CALL_SUBTEST( verify_component_wise(ref, ref); ); - CALL_SUBTEST( res = x.zeta(q); verify_component_wise(res, ref); ); - CALL_SUBTEST( res = zeta(x,q); verify_component_wise(res, ref); ); - } - - // digamma - { - ArrayType x(7), res(7), ref(7); - x << 1, 1.5, 4, -10.5, 10000.5, 0, -1; - ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, plusinf, plusinf; - CALL_SUBTEST( verify_component_wise(ref, ref); ); - - CALL_SUBTEST( res = x.digamma(); verify_component_wise(res, ref); ); - CALL_SUBTEST( res = digamma(x); verify_component_wise(res, ref); ); - } - - -#if EIGEN_HAS_C99_MATH - { - ArrayType n(11), x(11), res(11), ref(11); - n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170; - x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64; - ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927; - CALL_SUBTEST( verify_component_wise(ref, ref); ); - - if(sizeof(RealScalar)>=8) { // double - // Reason for commented line: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 - // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res, ref); ); - CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res, ref); ); - } - else { - // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res.head(8), ref.head(8)); ); - CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res.head(8), ref.head(8)); ); - } - } -#endif - -#if EIGEN_HAS_C99_MATH - { - // Inputs and ground truth generated with scipy via: - // a = np.logspace(-3, 3, 5) - 1e-3 - // b = np.logspace(-3, 3, 5) - 1e-3 - // x = np.linspace(-0.1, 1.1, 5) - // (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x)) - // full_a = full_a.flatten().tolist() # same for full_b, full_x - // v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist() - // - // Note in Eigen, we call betainc with arguments in the order (x, a, b). - ArrayType a(125); - ArrayType b(125); - ArrayType x(125); - ArrayType v(125); - ArrayType res(125); - - a << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, - 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, - 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, - 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, - 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, - 999.999, 999.999, 999.999; - - b << 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, - 0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, - 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, - 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, - 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, - 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, - 31.62177660168379, 31.62177660168379, 31.62177660168379, - 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, - 999.999, 999.999; - - x << -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, - 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, - 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, - 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, - -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, - 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, - 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, - 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, - 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, - 0.8, 1.1; - - v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, - nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, - nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan, - 0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan, - 0.999995949033062, 0.9999999999993698, 0.9999999999999999, nan, nan, - 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, - nan, nan, nan, nan, nan, 0.006827081192655869, 0.0210336989586256, - 0.04813160422599567, nan, nan, 0.20014344256217678, 0.5000000000000001, - 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403, - 0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999, - 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan, - 1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06, - nan, nan, 7.864342668429763e-23, 3.015969667594166e-10, - 0.0008598571564165444, nan, nan, 6.031987710123844e-08, - 0.5000000000000007, 0.9999999396801229, nan, nan, 0.9999999999999999, - 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan, - nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan, - 0.0, 9.275871147869727e-302, 1.2232913026152827e-97, nan, nan, 0.0, - 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan, - 2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan; - - CALL_SUBTEST(res = betainc(a, b, x); - verify_component_wise(res, v);); - } - - // Test various properties of betainc - { - ArrayType m1 = ArrayType::Random(32); - ArrayType m2 = ArrayType::Random(32); - ArrayType m3 = ArrayType::Random(32); - ArrayType one = ArrayType::Constant(32, Scalar(1.0)); - const Scalar eps = std::numeric_limits::epsilon(); - ArrayType a = (m1 * 4.0).exp(); - ArrayType b = (m2 * 4.0).exp(); - ArrayType x = m3.abs(); - - // betainc(a, 1, x) == x**a - CALL_SUBTEST( - ArrayType test = betainc(a, one, x); - ArrayType expected = x.pow(a); - verify_component_wise(test, expected);); - - // betainc(1, b, x) == 1 - (1 - x)**b - CALL_SUBTEST( - ArrayType test = betainc(one, b, x); - ArrayType expected = one - (one - x).pow(b); - verify_component_wise(test, expected);); - - // betainc(a, b, x) == 1 - betainc(b, a, 1-x) - CALL_SUBTEST( - ArrayType test = betainc(a, b, x) + betainc(b, a, one - x); - ArrayType expected = one; - verify_component_wise(test, expected);); - - // betainc(a+1, b, x) = betainc(a, b, x) - x**a * (1 - x)**b / (a * beta(a, b)) - CALL_SUBTEST( - ArrayType num = x.pow(a) * (one - x).pow(b); - ArrayType denom = a * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); - // Add eps to rhs and lhs so that component-wise test doesn't result in - // nans when both outputs are zeros. - ArrayType expected = betainc(a, b, x) - num / denom + eps; - ArrayType test = betainc(a + one, b, x) + eps; - if (sizeof(Scalar) >= 8) { // double - verify_component_wise(test, expected); - } else { - // Reason for limited test: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 - verify_component_wise(test.head(8), expected.head(8)); - }); - - // betainc(a, b+1, x) = betainc(a, b, x) + x**a * (1 - x)**b / (b * beta(a, b)) - CALL_SUBTEST( - // Add eps to rhs and lhs so that component-wise test doesn't result in - // nans when both outputs are zeros. - ArrayType num = x.pow(a) * (one - x).pow(b); - ArrayType denom = b * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); - ArrayType expected = betainc(a, b, x) + num / denom + eps; - ArrayType test = betainc(a, b + one, x) + eps; - verify_component_wise(test, expected);); - } -#endif -} - void test_array() { -#ifndef EIGEN_HAS_C99_MATH - std::cerr << "WARNING: testing of special math functions disabled" << std::endl; -#endif - for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( array(Array()) ); CALL_SUBTEST_2( array(Array22f()) ); @@ -812,7 +491,4 @@ void test_array() VERIFY((internal::is_same< internal::global_math_functions_filtering_base::type, ArrayBase >::value)); - - CALL_SUBTEST_7(array_special_functions()); - CALL_SUBTEST_7(array_special_functions()); } diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt index 6d0cf4f9d..7478b6b0d 100644 --- a/unsupported/Eigen/CMakeLists.txt +++ b/unsupported/Eigen/CMakeLists.txt @@ -17,6 +17,7 @@ set(Eigen_HEADERS Polynomials Skyline SparseExtra + SpecialFunctions Splines ) diff --git a/unsupported/Eigen/CXX11/Tensor b/unsupported/Eigen/CXX11/Tensor index 79bac2f67..f7b94cee1 100644 --- a/unsupported/Eigen/CXX11/Tensor +++ b/unsupported/Eigen/CXX11/Tensor @@ -15,6 +15,7 @@ #include +#include "../SpecialFunctions" #include "src/util/CXX11Meta.h" #include "src/util/MaxSizeVector.h" diff --git a/unsupported/Eigen/SpecialFunctions b/unsupported/Eigen/SpecialFunctions new file mode 100644 index 000000000..b4ed83b30 --- /dev/null +++ b/unsupported/Eigen/SpecialFunctions @@ -0,0 +1,57 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud +// +// 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_SPECIALFUNCTIONS_MODULE +#define EIGEN_SPECIALFUNCTIONS_MODULE + +#include "../../Eigen/Core" + +#include "../../Eigen/src/Core/util/DisableStupidWarnings.h" + +namespace Eigen { + +/** + * \defgroup SpecialFunctions_Module Special math functions module + * + * This module features additional coefficient-wise math functions available + * within the numext:: namespace for the scalar version, and as method and/or free + * functions of Array. Those include: + * + * - erf + * - erfc + * - lgamma + * - igamma + * - igammac + * - digamma + * - polygamma + * - zeta + * - betainc + * + * \code + * #include + * \endcode + */ +//@{ + +} + +#include "src/SpecialFunctions/SpecialFunctionsImpl.h" +#include "src/SpecialFunctions/SpecialFunctionsPacketMath.h" +#include "src/SpecialFunctions/SpecialFunctionsHalf.h" +#include "src/SpecialFunctions/SpecialFunctionsFunctors.h" +#include "src/SpecialFunctions/SpecialFunctionsArrayAPI.h" + +namespace Eigen { +//@} +} + + +#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_SPECIALFUNCTIONS_MODULE diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index a7e8c7553..f42946793 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -11,5 +11,6 @@ ADD_SUBDIRECTORY(NumericalDiff) ADD_SUBDIRECTORY(Polynomials) ADD_SUBDIRECTORY(Skyline) ADD_SUBDIRECTORY(SparseExtra) +ADD_SUBDIRECTORY(SpecialFunctions) ADD_SUBDIRECTORY(KroneckerProduct) ADD_SUBDIRECTORY(Splines) diff --git a/unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt b/unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt new file mode 100644 index 000000000..00c7afd7b --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_SpecialFunctions_SRCS "*.h") + +INSTALL(FILES + ${Eigen_SpecialFunctions_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/SpecialFunctions COMPONENT Devel + ) diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h new file mode 100644 index 000000000..ed415db99 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h @@ -0,0 +1,124 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud +// +// 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_SPECIALFUNCTIONS_ARRAYAPI_H +#define EIGEN_SPECIALFUNCTIONS_ARRAYAPI_H + +namespace Eigen { + +/** \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 +inline const Eigen::CwiseBinaryOp, const Derived, const ExponentDerived> +igamma(const Eigen::ArrayBase& a, const Eigen::ArrayBase& x) +{ + return Eigen::CwiseBinaryOp, 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 +inline const Eigen::CwiseBinaryOp, const Derived, const ExponentDerived> +igammac(const Eigen::ArrayBase& a, const Eigen::ArrayBase& x) +{ + return Eigen::CwiseBinaryOp, 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 +inline const Eigen::CwiseBinaryOp, const DerivedN, const DerivedX> +polygamma(const Eigen::ArrayBase& n, const Eigen::ArrayBase& x) +{ + return Eigen::CwiseBinaryOp, 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 +inline const Eigen::CwiseTernaryOp, const ArgADerived, const ArgBDerived, const ArgXDerived> +betainc(const Eigen::ArrayBase& a, const Eigen::ArrayBase& b, const Eigen::ArrayBase& x) +{ + return Eigen::CwiseTernaryOp, 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 +inline const Eigen::CwiseBinaryOp, const DerivedX, const DerivedQ> +zeta(const Eigen::ArrayBase& x, const Eigen::ArrayBase& q) +{ + return Eigen::CwiseBinaryOp, const DerivedX, const DerivedQ>( + x.derived(), + q.derived() + ); +} + +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_ARRAYAPI_H diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h new file mode 100644 index 000000000..d8f2363be --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h @@ -0,0 +1,236 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Eugene Brevdo +// Copyright (C) 2016 Gael Guennebaud +// +// 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_SPECIALFUNCTIONS_FUNCTORS_H +#define EIGEN_SPECIALFUNCTIONS_FUNCTORS_H + +namespace Eigen { + +namespace internal { + + +/** \internal + * \brief Template functor to compute the incomplete gamma function igamma(a, x) + * + * \sa class CwiseBinaryOp, Cwise::igamma + */ +template struct scalar_igamma_op : binary_op_base +{ + 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 + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const { + return internal::pigamma(a, x); + } +}; +template +struct functor_traits > { + enum { + // Guesstimate + Cost = 20 * NumTraits::MulCost + 10 * NumTraits::AddCost, + PacketAccess = packet_traits::HasIGamma + }; +}; + + +/** \internal + * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x) + * + * \sa class CwiseBinaryOp, Cwise::igammac + */ +template struct scalar_igammac_op : binary_op_base +{ + 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 + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const + { + return internal::pigammac(a, x); + } +}; +template +struct functor_traits > { + enum { + // Guesstimate + Cost = 20 * NumTraits::MulCost + 10 * NumTraits::AddCost, + PacketAccess = packet_traits::HasIGammac + }; +}; + + +/** \internal + * \brief Template functor to compute the incomplete beta integral betainc(a, b, x) + * + */ +template 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 + 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 +struct functor_traits > { + enum { + // Guesstimate + Cost = 400 * NumTraits::MulCost + 400 * NumTraits::AddCost, + PacketAccess = packet_traits::HasBetaInc + }; +}; + + +/** \internal + * \brief Template functor to compute the natural log of the absolute + * value of Gamma of a scalar + * \sa class CwiseUnaryOp, Cwise::lgamma() + */ +template 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::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::plgamma(a); } +}; +template +struct functor_traits > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, + PacketAccess = packet_traits::HasLGamma + }; +}; + +/** \internal + * \brief Template functor to compute psi, the derivative of lgamma of a scalar. + * \sa class CwiseUnaryOp, Cwise::digamma() + */ +template 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::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::pdigamma(a); } +}; +template +struct functor_traits > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, + PacketAccess = packet_traits::HasDiGamma + }; +}; + +/** \internal + * \brief Template functor to compute the Riemann Zeta function of two arguments. + * \sa class CwiseUnaryOp, Cwise::zeta() + */ +template 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::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); } +}; +template +struct functor_traits > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, + PacketAccess = packet_traits::HasZeta + }; +}; + +/** \internal + * \brief Template functor to compute the polygamma function. + * \sa class CwiseUnaryOp, Cwise::polygamma() + */ +template 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::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& n, const Packet& x) const { return internal::ppolygamma(n, x); } +}; +template +struct functor_traits > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, + PacketAccess = packet_traits::HasPolygamma + }; +}; + +/** \internal + * \brief Template functor to compute the Gauss error function of a + * scalar + * \sa class CwiseUnaryOp, Cwise::erf() + */ +template 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::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perf(a); } +}; +template +struct functor_traits > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, + PacketAccess = packet_traits::HasErf + }; +}; + +/** \internal + * \brief Template functor to compute the Complementary Error Function + * of a scalar + * \sa class CwiseUnaryOp, Cwise::erfc() + */ +template 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::type Packet; + EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::perfc(a); } +}; +template +struct functor_traits > +{ + enum { + // Guesstimate + Cost = 10 * NumTraits::MulCost + 5 * NumTraits::AddCost, + PacketAccess = packet_traits::HasErfc + }; +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_FUNCTORS_H diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h new file mode 100644 index 000000000..553bcda6a --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h @@ -0,0 +1,47 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// 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_SPECIALFUNCTIONS_HALF_H +#define EIGEN_SPECIALFUNCTIONS_HALF_H + +namespace Eigen { +namespace numext { + +#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(a))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half digamma(const Eigen::half& a) { + return Eigen::half(Eigen::numext::digamma(static_cast(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(x), static_cast(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(n), static_cast(x))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erf(const Eigen::half& a) { + return Eigen::half(Eigen::numext::erf(static_cast(a))); +} +template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::half& a) { + return Eigen::half(Eigen::numext::erfc(static_cast(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(a), static_cast(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(a), static_cast(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(a), static_cast(b), static_cast(x))); +} +#endif + +} // end namespace numext +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_HALF_H diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h new file mode 100644 index 000000000..52619fc0c --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -0,0 +1,1551 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2015 Eugene Brevdo +// +// 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( 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 +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::run(x, coef) * x + coef[N]; + } +}; + +template +struct polevl { + 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 +struct lgamma_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +template +struct lgamma_retval { + typedef Scalar type; +}; + +#if EIGEN_HAS_C99_MATH +template <> +struct lgamma_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE float run(float x) { return ::lgammaf(x); } +}; + +template <> +struct lgamma_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE double run(double x) { return ::lgamma(x); } +}; +#endif + +/**************************************************************************** + * Implementation of digamma (psi), based on Cephes * + ****************************************************************************/ + +template +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 +struct digamma_impl_maybe_poly { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + + +template <> +struct digamma_impl_maybe_poly { + 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::run(z, A); + } else return 0.0f; + } +}; + +template <> +struct digamma_impl_maybe_poly { + 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::run(z, A); + } + else return 0.0; + } +}; + +template +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::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::run(s); + + y = numext::log(s) - (half / s) - y - w; + + return (negative) ? y - nz : y; + } +}; + +/**************************************************************************** + * Implementation of erf, requires C++11/C99 * + ****************************************************************************/ + +template +struct erf_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +template +struct erf_retval { + typedef Scalar type; +}; + +#if EIGEN_HAS_C99_MATH +template <> +struct erf_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); } +}; + +template <> +struct erf_impl { + 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 +struct erfc_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +template +struct erfc_retval { + typedef Scalar type; +}; + +#if EIGEN_HAS_C99_MATH +template <> +struct erfc_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); } +}; + +template <> +struct erfc_impl { + 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 +struct igammac_retval { + typedef Scalar type; +}; + +// NOTE: cephes_helper is also used to implement zeta +template +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 { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE float machep() { + return NumTraits::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::epsilon() / 2); + } + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE float biginv() { + // epsneg + return machep(); + } +}; + +template <> +struct cephes_helper { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE double machep() { + return NumTraits::epsilon() / 2; // 1.0 - machep == 1.0 + } + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE double big() { + return 1.0 / NumTraits::epsilon(); + } + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE double biginv() { + // inverse of eps + return NumTraits::epsilon(); + } +}; + +#if !EIGEN_HAS_C99_MATH + +template +struct igammac_impl { + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +#else + +template struct igamma_impl; // predeclare igamma_impl + +template +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::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::Impl(a, x)); + } + + return Impl(a, x); + } + + private: + /* igamma_impl calls igammac_impl::Impl. */ + friend struct igamma_impl; + + /* 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::machep(); + const Scalar maxlog = numext::log(NumTraits::highest()); + const Scalar big = cephes_helper::big(); + const Scalar biginv = cephes_helper::biginv(); + const Scalar inf = NumTraits::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::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 +struct igamma_retval { + typedef Scalar type; +}; + +#if !EIGEN_HAS_C99_MATH + +template +struct igamma_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +#else + +template +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::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::Impl(a, x)); + } + + return Impl(a, x); + } + + private: + /* igammac_impl calls igamma_impl::Impl. */ + friend struct igammac_impl; + + /* 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::machep(); + const Scalar maxlog = numext::log(NumTraits::highest()); + + Scalar ans, ax, c, r; + + /* Compute x**a * exp(-x) / gamma(a) */ + ax = a * numext::log(x) - x - lgamma_impl::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 +struct zeta_retval { + typedef Scalar type; +}; + +template +struct zeta_impl_series { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +template <> +struct zeta_impl_series { + 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 { + 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 +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::infinity(); + const Scalar zero = 0.0, half = 0.5, one = 1.0; + const Scalar machep = cephes_helper::machep(); + const Scalar nan = NumTraits::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::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 +struct polygamma_retval { + typedef Scalar type; +}; + +#if !EIGEN_HAS_C99_MATH + +template +struct polygamma_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +#else + +template +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::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::run(x); + } + // Use the same implementation as scipy + else { + Scalar factorial = numext::exp(lgamma_impl::run(nplus)); + return numext::pow(-one, nplus) * factorial * zeta_impl::run(nplus, x); + } + } +}; + +#endif // EIGEN_HAS_C99_MATH + +/************************************************************************************************ + * Implementation of betainc (incomplete beta integral), based on Cephes but requires C++11/C99 * + ************************************************************************************************/ + +template +struct betainc_retval { + typedef Scalar type; +}; + +#if !EIGEN_HAS_C99_MATH + +template +struct betainc_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) { + EIGEN_STATIC_ASSERT((internal::is_same::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +#else + +template +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::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 +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::value || + internal::is_same::value), + THIS_TYPE_IS_NOT_SUPPORTED); + const Scalar big = cephes_helper::big(); + const Scalar machep = cephes_helper::machep(); + const Scalar biginv = cephes_helper::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::value) ? 100 : 300; + const Scalar thresh = + (internal::is_same::value) ? machep : Scalar(3) * machep; + Scalar r = (internal::is_same::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 +struct betainc_helper {}; + +template <> +struct betainc_helper { + /* 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::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::run(a, b, x, true /* small_branch */); + t = b * numext::log(t); + } else { + ans = incbeta_cfe::run(a, b, x, false /* small_branch */); + t = (b - 1.0f) * numext::log(t); + } + + t += a * numext::log(x) + lgamma_impl::run(a + b) - + lgamma_impl::run(a) - lgamma_impl::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::machep(); + + y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a); + y -= lgamma_impl::run(a) + lgamma_impl::run(b); + y += lgamma_impl::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 { + EIGEN_DEVICE_FUNC + static float run(float a, float b, float x) { + const float nan = NumTraits::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::incbsa(a + 1.0f, b, x); + t = a * numext::log(x) + b * numext::log1p(-x) + + lgamma_impl::run(a + b) - lgamma_impl::run(a + 1.0f) - + lgamma_impl::run(b); + return (ans + numext::exp(t)); + } else { + return betainc_helper::incbsa(a, b, x); + } + } +}; + +template <> +struct betainc_helper { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE double incbps(double a, double b, double x) { + const double machep = cephes_helper::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::run(a + b) - lgamma_impl::run(a) - + lgamma_impl::run(b) + u + numext::log(s); + return s = numext::exp(t); + } +}; + +template <> +struct betainc_impl { + EIGEN_DEVICE_FUNC + static double run(double aa, double bb, double xx) { + const double nan = NumTraits::quiet_NaN(); + const double machep = cephes_helper::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::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::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::run(a, b, x, true /* small_branch */); + } else { + w = incbeta_cfe::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::run(a + b) - lgamma_impl::run(a) - + lgamma_impl::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 +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar) + lgamma(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar) + digamma(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x); +} + +template +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 +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 +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar) + erf(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x); +} + +template +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) + erfc(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x); +} + +template +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 +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 +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/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h new file mode 100644 index 000000000..46d60d323 --- /dev/null +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h @@ -0,0 +1,58 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud +// +// 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_SPECIALFUNCTIONS_PACKETMATH_H +#define EIGEN_SPECIALFUNCTIONS_PACKETMATH_H + +namespace Eigen { + +namespace internal { + +/** \internal \returns the ln(|gamma(\a a)|) (coeff-wise) */ +template 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 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 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 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 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 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 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 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 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); } + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SPECIALFUNCTIONS_PACKETMATH_H + diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index c9a70d7a7..5137b51cf 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -109,6 +109,7 @@ ei_add_test(gmres) ei_add_test(minres) ei_add_test(levenberg_marquardt) ei_add_test(kronecker_product) +ei_add_test(special_functions) # TODO: The following test names are prefixed with the cxx11 string, since historically # the tests depended on c++11. This isn't the case anymore so we ought to rename them. diff --git a/unsupported/test/special_functions.cpp b/unsupported/test/special_functions.cpp new file mode 100644 index 000000000..057fb3e92 --- /dev/null +++ b/unsupported/test/special_functions.cpp @@ -0,0 +1,345 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include "../Eigen/SpecialFunctions" + +template +void verify_component_wise(const X& x, const Y& y) +{ + for(Index i=0; i void array_special_functions() +{ + using std::abs; + using std::sqrt; + typedef typename ArrayType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + Scalar plusinf = std::numeric_limits::infinity(); + Scalar nan = std::numeric_limits::quiet_NaN(); + + Index rows = internal::random(1,30); + Index cols = 1; + + // API + { + ArrayType m1 = ArrayType::Random(rows,cols); +#if EIGEN_HAS_C99_MATH + VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1)); + VERIFY_IS_APPROX(m1.digamma(), digamma(m1)); + VERIFY_IS_APPROX(m1.erf(), erf(m1)); + VERIFY_IS_APPROX(m1.erfc(), erfc(m1)); +#endif // EIGEN_HAS_C99_MATH + } + + +#if EIGEN_HAS_C99_MATH + // check special functions (comparing against numpy implementation) + if (!NumTraits::IsComplex) + { + + { + ArrayType m1 = ArrayType::Random(rows,cols); + ArrayType m2 = ArrayType::Random(rows,cols); + + // Test various propreties of igamma & igammac. These are normalized + // gamma integrals where + // igammac(a, x) = Gamma(a, x) / Gamma(a) + // igamma(a, x) = gamma(a, x) / Gamma(a) + // where Gamma and gamma are considered the standard unnormalized + // upper and lower incomplete gamma functions, respectively. + ArrayType a = m1.abs() + 2; + ArrayType x = m2.abs() + 2; + ArrayType zero = ArrayType::Zero(rows, cols); + ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0)); + ArrayType a_m1 = a - one; + ArrayType Gamma_a_x = Eigen::igammac(a, x) * a.lgamma().exp(); + ArrayType Gamma_a_m1_x = Eigen::igammac(a_m1, x) * a_m1.lgamma().exp(); + ArrayType gamma_a_x = Eigen::igamma(a, x) * a.lgamma().exp(); + ArrayType gamma_a_m1_x = Eigen::igamma(a_m1, x) * a_m1.lgamma().exp(); + + // Gamma(a, 0) == Gamma(a) + VERIFY_IS_APPROX(Eigen::igammac(a, zero), one); + + // Gamma(a, x) + gamma(a, x) == Gamma(a) + VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp()); + + // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) + VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).exp()); + + // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) + VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).exp()); + } + + { + // Check exact values of igamma and igammac against a third party calculation. + Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; + Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; + + // location i*6+j corresponds to a_s[i], x_s[j]. + Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan}, + {0.0, 0.6321205588285578, 0.7768698398515702, + 0.9816843611112658, 9.999500016666262e-05, 1.0}, + {0.0, 0.4275932955291202, 0.608374823728911, + 0.9539882943107686, 7.522076445089201e-07, 1.0}, + {0.0, 0.01898815687615381, 0.06564245437845008, + 0.5665298796332909, 4.166333347221828e-18, 1.0}, + {0.0, 0.9999780593618628, 0.9999899967080838, + 0.9999996219837988, 0.9991370418689945, 1.0}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}}; + Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan}, + {1.0, 0.36787944117144233, 0.22313016014842982, + 0.018315638888734182, 0.9999000049998333, 0.0}, + {1.0, 0.5724067044708798, 0.3916251762710878, + 0.04601170568923136, 0.9999992477923555, 0.0}, + {1.0, 0.9810118431238462, 0.9343575456215499, + 0.4334701203667089, 1.0, 0.0}, + {1.0, 2.1940638138146658e-05, 1.0003291916285e-05, + 3.7801620118431334e-07, 0.0008629581310054535, + 0.0}, + {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}}; + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 6; ++j) { + if ((std::isnan)(igamma_s[i][j])) { + VERIFY((std::isnan)(numext::igamma(a_s[i], x_s[j]))); + } else { + VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]); + } + + if ((std::isnan)(igammac_s[i][j])) { + VERIFY((std::isnan)(numext::igammac(a_s[i], x_s[j]))); + } else { + VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]); + } + } + } + } + } +#endif // EIGEN_HAS_C99_MATH + + // Check the zeta function against scipy.special.zeta + { + ArrayType x(7), q(7), res(7), ref(7); + x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9; + q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345; + ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan; + CALL_SUBTEST( verify_component_wise(ref, ref); ); + CALL_SUBTEST( res = x.zeta(q); verify_component_wise(res, ref); ); + CALL_SUBTEST( res = zeta(x,q); verify_component_wise(res, ref); ); + } + + // digamma + { + ArrayType x(7), res(7), ref(7); + x << 1, 1.5, 4, -10.5, 10000.5, 0, -1; + ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, plusinf, plusinf; + CALL_SUBTEST( verify_component_wise(ref, ref); ); + + CALL_SUBTEST( res = x.digamma(); verify_component_wise(res, ref); ); + CALL_SUBTEST( res = digamma(x); verify_component_wise(res, ref); ); + } + + +#if EIGEN_HAS_C99_MATH + { + ArrayType n(11), x(11), res(11), ref(11); + n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170; + x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64; + ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927; + CALL_SUBTEST( verify_component_wise(ref, ref); ); + + if(sizeof(RealScalar)>=8) { // double + // Reason for commented line: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 + // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res, ref); ); + CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res, ref); ); + } + else { + // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res.head(8), ref.head(8)); ); + CALL_SUBTEST( res = polygamma(n,x); verify_component_wise(res.head(8), ref.head(8)); ); + } + } +#endif + +#if EIGEN_HAS_C99_MATH + { + // Inputs and ground truth generated with scipy via: + // a = np.logspace(-3, 3, 5) - 1e-3 + // b = np.logspace(-3, 3, 5) - 1e-3 + // x = np.linspace(-0.1, 1.1, 5) + // (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x)) + // full_a = full_a.flatten().tolist() # same for full_b, full_x + // v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist() + // + // Note in Eigen, we call betainc with arguments in the order (x, a, b). + ArrayType a(125); + ArrayType b(125); + ArrayType x(125); + ArrayType v(125); + ArrayType res(125); + + a << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, + 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, + 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, + 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, + 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, + 999.999, 999.999, 999.999; + + b << 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, + 0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, + 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, + 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, + 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, + 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, + 31.62177660168379, 31.62177660168379, 31.62177660168379, + 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, + 999.999, 999.999; + + x << -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, + 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, + 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, + 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, + -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, + 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, + 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, + 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, + 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, + 0.8, 1.1; + + v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, + nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, + nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan, + 0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan, + 0.999995949033062, 0.9999999999993698, 0.9999999999999999, nan, nan, + 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, + nan, nan, nan, nan, nan, 0.006827081192655869, 0.0210336989586256, + 0.04813160422599567, nan, nan, 0.20014344256217678, 0.5000000000000001, + 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403, + 0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999, + 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan, + 1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06, + nan, nan, 7.864342668429763e-23, 3.015969667594166e-10, + 0.0008598571564165444, nan, nan, 6.031987710123844e-08, + 0.5000000000000007, 0.9999999396801229, nan, nan, 0.9999999999999999, + 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan, + nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan, + 0.0, 9.275871147869727e-302, 1.2232913026152827e-97, nan, nan, 0.0, + 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan, + 2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan; + + CALL_SUBTEST(res = betainc(a, b, x); + verify_component_wise(res, v);); + } + + // Test various properties of betainc + { + ArrayType m1 = ArrayType::Random(32); + ArrayType m2 = ArrayType::Random(32); + ArrayType m3 = ArrayType::Random(32); + ArrayType one = ArrayType::Constant(32, Scalar(1.0)); + const Scalar eps = std::numeric_limits::epsilon(); + ArrayType a = (m1 * 4.0).exp(); + ArrayType b = (m2 * 4.0).exp(); + ArrayType x = m3.abs(); + + // betainc(a, 1, x) == x**a + CALL_SUBTEST( + ArrayType test = betainc(a, one, x); + ArrayType expected = x.pow(a); + verify_component_wise(test, expected);); + + // betainc(1, b, x) == 1 - (1 - x)**b + CALL_SUBTEST( + ArrayType test = betainc(one, b, x); + ArrayType expected = one - (one - x).pow(b); + verify_component_wise(test, expected);); + + // betainc(a, b, x) == 1 - betainc(b, a, 1-x) + CALL_SUBTEST( + ArrayType test = betainc(a, b, x) + betainc(b, a, one - x); + ArrayType expected = one; + verify_component_wise(test, expected);); + + // betainc(a+1, b, x) = betainc(a, b, x) - x**a * (1 - x)**b / (a * beta(a, b)) + CALL_SUBTEST( + ArrayType num = x.pow(a) * (one - x).pow(b); + ArrayType denom = a * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); + // Add eps to rhs and lhs so that component-wise test doesn't result in + // nans when both outputs are zeros. + ArrayType expected = betainc(a, b, x) - num / denom + eps; + ArrayType test = betainc(a + one, b, x) + eps; + if (sizeof(Scalar) >= 8) { // double + verify_component_wise(test, expected); + } else { + // Reason for limited test: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 + verify_component_wise(test.head(8), expected.head(8)); + }); + + // betainc(a, b+1, x) = betainc(a, b, x) + x**a * (1 - x)**b / (b * beta(a, b)) + CALL_SUBTEST( + // Add eps to rhs and lhs so that component-wise test doesn't result in + // nans when both outputs are zeros. + ArrayType num = x.pow(a) * (one - x).pow(b); + ArrayType denom = b * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp(); + ArrayType expected = betainc(a, b, x) + num / denom + eps; + ArrayType test = betainc(a, b + one, x) + eps; + verify_component_wise(test, expected);); + } +#endif +} + +void test_special_functions() +{ + CALL_SUBTEST_1(array_special_functions()); + CALL_SUBTEST_2(array_special_functions()); +} -- cgit v1.2.3