diff options
Diffstat (limited to 'unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h')
| -rw-r--r-- | unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h | 1048 |
1 files changed, 764 insertions, 284 deletions
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index f524d7137..f1c260e29 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -36,66 +36,6 @@ namespace internal { // Good luck with your project, // Steve -namespace cephes { - -/* polevl (modified for Eigen) - * - * Evaluate polynomial - * - * - * - * SYNOPSIS: - * - * int N; - * Scalar x, y, coef[N+1]; - * - * y = polevl<decltype(x), N>( x, coef); - * - * - * - * DESCRIPTION: - * - * Evaluates polynomial of degree N: - * - * 2 N - * y = C + C x + C x +...+ C x - * 0 1 2 N - * - * Coefficients are stored in reverse order: - * - * coef[0] = C , ..., coef[N] = C . - * N 0 - * - * The function p1evl() assumes that coef[N] = 1.0 and is - * omitted from the array. Its calling arguments are - * otherwise the same as polevl(). - * - * - * The Eigen implementation is templatized. For best speed, store - * coef as a const array (constexpr), e.g. - * - * const double coef[] = {1.0, 2.0, 3.0, ...}; - * - */ -template <typename Scalar, int N> -struct polevl { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar x, const Scalar coef[]) { - EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE); - - return polevl<Scalar, N - 1>::run(x, coef) * x + coef[N]; - } -}; - -template <typename Scalar> -struct polevl<Scalar, 0> { - EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar, const Scalar coef[]) { - return coef[0]; - } -}; - -} // end namespace cephes /**************************************************************************** * Implementation of lgamma, requires C++11/C99 * @@ -117,13 +57,27 @@ struct lgamma_retval { }; #if EIGEN_HAS_C99_MATH +// Since glibc 2.19 +#if defined(__GLIBC__) && ((__GLIBC__>=2 && __GLIBC_MINOR__ >= 19) || __GLIBC__>2) \ + && (defined(_DEFAULT_SOURCE) || defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) +#define EIGEN_HAS_LGAMMA_R +#endif + +// Glibc versions before 2.19 +#if defined(__GLIBC__) && ((__GLIBC__==2 && __GLIBC_MINOR__ < 19) || __GLIBC__<2) \ + && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) +#define EIGEN_HAS_LGAMMA_R +#endif + template <> struct lgamma_impl<float> { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE float run(float x) { -#if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) - int signgam; - return ::lgammaf_r(x, &signgam); +#if !defined(EIGEN_GPU_COMPILE_PHASE) && defined (EIGEN_HAS_LGAMMA_R) && !defined(__APPLE__) + int dummy; + return ::lgammaf_r(x, &dummy); +#elif defined(SYCL_DEVICE_ONLY) + return cl::sycl::lgamma(x); #else return ::lgammaf(x); #endif @@ -134,14 +88,18 @@ template <> struct lgamma_impl<double> { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE double run(double x) { -#if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) - int signgam; - return ::lgamma_r(x, &signgam); +#if !defined(EIGEN_GPU_COMPILE_PHASE) && defined(EIGEN_HAS_LGAMMA_R) && !defined(__APPLE__) + int dummy; + return ::lgamma_r(x, &dummy); +#elif defined(SYCL_DEVICE_ONLY) + return cl::sycl::lgamma(x); #else return ::lgamma(x); #endif } }; + +#undef EIGEN_HAS_LGAMMA_R #endif /**************************************************************************** @@ -191,7 +149,7 @@ struct digamma_impl_maybe_poly<float> { float z; if (s < 1.0e8f) { z = 1.0f / (s * s); - return z * cephes::polevl<float, 3>::run(z, A); + return z * internal::ppolevl<float, 3>::run(z, A); } else return 0.0f; } }; @@ -213,7 +171,7 @@ struct digamma_impl_maybe_poly<double> { double z; if (s < 1.0e17) { z = 1.0 / (s * s); - return z * cephes::polevl<double, 6>::run(z, A); + return z * internal::ppolevl<double, 6>::run(z, A); } else return 0.0; } @@ -283,7 +241,7 @@ struct digamma_impl { Scalar p, q, nz, s, w, y; bool negative = false; - const Scalar maxnum = NumTraits<Scalar>::infinity(); + const Scalar nan = NumTraits<Scalar>::quiet_NaN(); const Scalar m_pi = Scalar(EIGEN_PI); const Scalar zero = Scalar(0); @@ -296,7 +254,7 @@ struct digamma_impl { q = x; p = numext::floor(q); if (p == q) { - return maxnum; + return nan; } /* Remove the zeros of tan(m_pi x) * by subtracting the nearest integer from x @@ -335,13 +293,63 @@ struct digamma_impl { * Implementation of erf, requires C++11/C99 * ****************************************************************************/ -template <typename Scalar> +/** \internal \returns the error function of \a a (coeff-wise) + Doesn't do anything fancy, just a 13/8-degree rational interpolant which + is accurate up to a couple of ulp in the range [-4, 4], outside of which + fl(erf(x)) = +/-1. + + This implementation works on both scalars and Ts. +*/ +template <typename T> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erf_float(const T& a_x) { + // Clamp the inputs to the range [-4, 4] since anything outside + // this range is +/-1.0f in single-precision. + const T plus_4 = pset1<T>(4.f); + const T minus_4 = pset1<T>(-4.f); + const T x = pmax(pmin(a_x, plus_4), minus_4); + // The monomial coefficients of the numerator polynomial (odd). + const T alpha_1 = pset1<T>(-1.60960333262415e-02f); + const T alpha_3 = pset1<T>(-2.95459980854025e-03f); + const T alpha_5 = pset1<T>(-7.34990630326855e-04f); + const T alpha_7 = pset1<T>(-5.69250639462346e-05f); + const T alpha_9 = pset1<T>(-2.10102402082508e-06f); + const T alpha_11 = pset1<T>(2.77068142495902e-08f); + const T alpha_13 = pset1<T>(-2.72614225801306e-10f); + + // The monomial coefficients of the denominator polynomial (even). + const T beta_0 = pset1<T>(-1.42647390514189e-02f); + const T beta_2 = pset1<T>(-7.37332916720468e-03f); + const T beta_4 = pset1<T>(-1.68282697438203e-03f); + const T beta_6 = pset1<T>(-2.13374055278905e-04f); + const T beta_8 = pset1<T>(-1.45660718464996e-05f); + + // Since the polynomials are odd/even, we need x^2. + const T x2 = pmul(x, x); + + // Evaluate the numerator polynomial p. + T p = pmadd(x2, alpha_13, alpha_11); + p = pmadd(x2, p, alpha_9); + p = pmadd(x2, p, alpha_7); + p = pmadd(x2, p, alpha_5); + p = pmadd(x2, p, alpha_3); + p = pmadd(x2, p, alpha_1); + p = pmul(x, p); + + // Evaluate the denominator polynomial p. + T q = pmadd(x2, beta_8, beta_6); + q = pmadd(x2, q, beta_4); + q = pmadd(x2, q, beta_2); + q = pmadd(x2, q, beta_0); + + // Divide the numerator by the denominator. + return pdiv(p, q); +} + +template <typename T> struct erf_impl { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE Scalar run(const Scalar) { - EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), - THIS_TYPE_IS_NOT_SUPPORTED); - return Scalar(0); + static EIGEN_STRONG_INLINE T run(const T& x) { + return generic_fast_erf_float(x); } }; @@ -354,13 +362,25 @@ struct erf_retval { template <> struct erf_impl<float> { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); } + static EIGEN_STRONG_INLINE float run(float x) { +#if defined(SYCL_DEVICE_ONLY) + return cl::sycl::erf(x); +#else + return generic_fast_erf_float(x); +#endif + } }; template <> struct erf_impl<double> { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double run(double x) { return ::erf(x); } + static EIGEN_STRONG_INLINE double run(double x) { +#if defined(SYCL_DEVICE_ONLY) + return cl::sycl::erf(x); +#else + return ::erf(x); +#endif + } }; #endif // EIGEN_HAS_C99_MATH @@ -387,16 +407,270 @@ struct erfc_retval { template <> struct erfc_impl<float> { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); } + static EIGEN_STRONG_INLINE float run(const float x) { +#if defined(SYCL_DEVICE_ONLY) + return cl::sycl::erfc(x); +#else + return ::erfcf(x); +#endif + } }; template <> struct erfc_impl<double> { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE double run(const double x) { return ::erfc(x); } + static EIGEN_STRONG_INLINE double run(const double x) { +#if defined(SYCL_DEVICE_ONLY) + return cl::sycl::erfc(x); +#else + return ::erfc(x); +#endif + } +}; +#endif // EIGEN_HAS_C99_MATH + + +/*************************************************************************** +* Implementation of ndtri. * +****************************************************************************/ + +/* Inverse of Normal distribution function (modified for Eigen). + * + * + * SYNOPSIS: + * + * double x, y, ndtri(); + * + * x = ndtri( y ); + * + * + * + * DESCRIPTION: + * + * Returns the argument, x, for which the area under the + * Gaussian probability density function (integrated from + * minus infinity to x) is equal to y. + * + * + * For small arguments 0 < y < exp(-2), the program computes + * z = sqrt( -2.0 * log(y) ); then the approximation is + * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). + * There are two rational functions P/Q, one for 0 < y < exp(-32) + * and the other for y up to exp(-2). For larger arguments, + * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0.125, 1 5500 9.5e-17 2.1e-17 + * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17 + * IEEE 0.125, 1 20000 7.2e-16 1.3e-16 + * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * ndtri domain x <= 0 -MAXNUM + * ndtri domain x >= 1 MAXNUM + * + */ + /* + 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 + */ + + +// TODO: Add a cheaper approximation for float. + + +template<typename T> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T flipsign( + const T& should_flipsign, const T& x) { + typedef typename unpacket_traits<T>::type Scalar; + const T sign_mask = pset1<T>(Scalar(-0.0)); + T sign_bit = pand<T>(should_flipsign, sign_mask); + return pxor<T>(sign_bit, x); +} + +template<> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double flipsign<double>( + const double& should_flipsign, const double& x) { + return should_flipsign == 0 ? x : -x; +} + +template<> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float flipsign<float>( + const float& should_flipsign, const float& x) { + return should_flipsign == 0 ? x : -x; +} + +// We split this computation in to two so that in the scalar path +// only one branch is evaluated (due to our template specialization of pselect +// being an if statement.) + +template <typename T, typename ScalarType> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_ndtri_gt_exp_neg_two(const T& b) { + const ScalarType p0[] = { + ScalarType(-5.99633501014107895267e1), + ScalarType(9.80010754185999661536e1), + ScalarType(-5.66762857469070293439e1), + ScalarType(1.39312609387279679503e1), + ScalarType(-1.23916583867381258016e0) + }; + const ScalarType q0[] = { + ScalarType(1.0), + ScalarType(1.95448858338141759834e0), + ScalarType(4.67627912898881538453e0), + ScalarType(8.63602421390890590575e1), + ScalarType(-2.25462687854119370527e2), + ScalarType(2.00260212380060660359e2), + ScalarType(-8.20372256168333339912e1), + ScalarType(1.59056225126211695515e1), + ScalarType(-1.18331621121330003142e0) + }; + const T sqrt2pi = pset1<T>(ScalarType(2.50662827463100050242e0)); + const T half = pset1<T>(ScalarType(0.5)); + T c, c2, ndtri_gt_exp_neg_two; + + c = psub(b, half); + c2 = pmul(c, c); + ndtri_gt_exp_neg_two = pmadd(c, pmul( + c2, pdiv( + internal::ppolevl<T, 4>::run(c2, p0), + internal::ppolevl<T, 8>::run(c2, q0))), c); + return pmul(ndtri_gt_exp_neg_two, sqrt2pi); +} + +template <typename T, typename ScalarType> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_ndtri_lt_exp_neg_two( + const T& b, const T& should_flipsign) { + /* Approximation for interval z = sqrt(-2 log a ) between 2 and 8 + * i.e., a between exp(-2) = .135 and exp(-32) = 1.27e-14. + */ + const ScalarType p1[] = { + ScalarType(4.05544892305962419923e0), + ScalarType(3.15251094599893866154e1), + ScalarType(5.71628192246421288162e1), + ScalarType(4.40805073893200834700e1), + ScalarType(1.46849561928858024014e1), + ScalarType(2.18663306850790267539e0), + ScalarType(-1.40256079171354495875e-1), + ScalarType(-3.50424626827848203418e-2), + ScalarType(-8.57456785154685413611e-4) + }; + const ScalarType q1[] = { + ScalarType(1.0), + ScalarType(1.57799883256466749731e1), + ScalarType(4.53907635128879210584e1), + ScalarType(4.13172038254672030440e1), + ScalarType(1.50425385692907503408e1), + ScalarType(2.50464946208309415979e0), + ScalarType(-1.42182922854787788574e-1), + ScalarType(-3.80806407691578277194e-2), + ScalarType(-9.33259480895457427372e-4) + }; + /* Approximation for interval z = sqrt(-2 log a ) between 8 and 64 + * i.e., a between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890. + */ + const ScalarType p2[] = { + ScalarType(3.23774891776946035970e0), + ScalarType(6.91522889068984211695e0), + ScalarType(3.93881025292474443415e0), + ScalarType(1.33303460815807542389e0), + ScalarType(2.01485389549179081538e-1), + ScalarType(1.23716634817820021358e-2), + ScalarType(3.01581553508235416007e-4), + ScalarType(2.65806974686737550832e-6), + ScalarType(6.23974539184983293730e-9) + }; + const ScalarType q2[] = { + ScalarType(1.0), + ScalarType(6.02427039364742014255e0), + ScalarType(3.67983563856160859403e0), + ScalarType(1.37702099489081330271e0), + ScalarType(2.16236993594496635890e-1), + ScalarType(1.34204006088543189037e-2), + ScalarType(3.28014464682127739104e-4), + ScalarType(2.89247864745380683936e-6), + ScalarType(6.79019408009981274425e-9) + }; + const T eight = pset1<T>(ScalarType(8.0)); + const T one = pset1<T>(ScalarType(1)); + const T neg_two = pset1<T>(ScalarType(-2)); + T x, x0, x1, z; + + x = psqrt(pmul(neg_two, plog(b))); + x0 = psub(x, pdiv(plog(x), x)); + z = pdiv(one, x); + x1 = pmul( + z, pselect( + pcmp_lt(x, eight), + pdiv(internal::ppolevl<T, 8>::run(z, p1), + internal::ppolevl<T, 8>::run(z, q1)), + pdiv(internal::ppolevl<T, 8>::run(z, p2), + internal::ppolevl<T, 8>::run(z, q2)))); + return flipsign(should_flipsign, psub(x0, x1)); +} + +template <typename T, typename ScalarType> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE +T generic_ndtri(const T& a) { + const T maxnum = pset1<T>(NumTraits<ScalarType>::infinity()); + const T neg_maxnum = pset1<T>(-NumTraits<ScalarType>::infinity()); + + const T zero = pset1<T>(ScalarType(0)); + const T one = pset1<T>(ScalarType(1)); + // exp(-2) + const T exp_neg_two = pset1<T>(ScalarType(0.13533528323661269189)); + T b, ndtri, should_flipsign; + + should_flipsign = pcmp_le(a, psub(one, exp_neg_two)); + b = pselect(should_flipsign, a, psub(one, a)); + + ndtri = pselect( + pcmp_lt(exp_neg_two, b), + generic_ndtri_gt_exp_neg_two<T, ScalarType>(b), + generic_ndtri_lt_exp_neg_two<T, ScalarType>(b, should_flipsign)); + + return pselect( + pcmp_le(a, zero), neg_maxnum, + pselect(pcmp_le(one, a), maxnum, ndtri)); +} + +template <typename Scalar> +struct ndtri_retval { + typedef Scalar type; +}; + +#if !EIGEN_HAS_C99_MATH + +template <typename Scalar> +struct ndtri_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar) { + EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), + THIS_TYPE_IS_NOT_SUPPORTED); + return Scalar(0); + } +}; + +# else + +template <typename Scalar> +struct ndtri_impl { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE Scalar run(const Scalar x) { + return generic_ndtri<Scalar, Scalar>(x); + } }; + #endif // EIGEN_HAS_C99_MATH + /************************************************************************************************************** * Implementation of igammac (complemented incomplete gamma integral), based on Cephes but requires C++11/C99 * **************************************************************************************************************/ @@ -452,6 +726,228 @@ struct cephes_helper<double> { } }; +enum IgammaComputationMode { VALUE, DERIVATIVE, SAMPLE_DERIVATIVE }; + +template <typename Scalar> +EIGEN_DEVICE_FUNC +static EIGEN_STRONG_INLINE Scalar main_igamma_term(Scalar a, Scalar x) { + /* Compute x**a * exp(-x) / gamma(a) */ + Scalar logax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); + if (logax < -numext::log(NumTraits<Scalar>::highest()) || + // Assuming x and a aren't Nan. + (numext::isnan)(logax)) { + return Scalar(0); + } + return numext::exp(logax); +} + +template <typename Scalar, IgammaComputationMode mode> +EIGEN_DEVICE_FUNC +int igamma_num_iterations() { + /* Returns the maximum number of internal iterations for igamma computation. + */ + if (mode == VALUE) { + return 2000; + } + + if (internal::is_same<Scalar, float>::value) { + return 200; + } else if (internal::is_same<Scalar, double>::value) { + return 500; + } else { + return 2000; + } +} + +template <typename Scalar, IgammaComputationMode mode> +struct igammac_cf_impl { + /* Computes igamc(a, x) or derivative (depending on the mode) + * using the continued fraction expansion of the complementary + * incomplete Gamma function. + * + * Preconditions: + * a > 0 + * x >= 1 + * x >= a + */ + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + const Scalar zero = 0; + const Scalar one = 1; + const Scalar two = 2; + const Scalar machep = cephes_helper<Scalar>::machep(); + const Scalar big = cephes_helper<Scalar>::big(); + const Scalar biginv = cephes_helper<Scalar>::biginv(); + + if ((numext::isinf)(x)) { + return zero; + } + + Scalar ax = main_igamma_term<Scalar>(a, x); + // This is independent of mode. If this value is zero, + // then the function value is zero. If the function value is zero, + // then we are in a neighborhood where the function value evalutes to zero, + // so the derivative is zero. + if (ax == zero) { + return zero; + } + + // continued fraction + Scalar y = one - a; + Scalar z = x + y + one; + Scalar c = zero; + Scalar pkm2 = one; + Scalar qkm2 = x; + Scalar pkm1 = x + one; + Scalar qkm1 = z * x; + Scalar ans = pkm1 / qkm1; + + Scalar dpkm2_da = zero; + Scalar dqkm2_da = zero; + Scalar dpkm1_da = zero; + Scalar dqkm1_da = -x; + Scalar dans_da = (dpkm1_da - ans * dqkm1_da) / qkm1; + + for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) { + c += one; + y += one; + z += two; + + Scalar yc = y * c; + Scalar pk = pkm1 * z - pkm2 * yc; + Scalar qk = qkm1 * z - qkm2 * yc; + + Scalar dpk_da = dpkm1_da * z - pkm1 - dpkm2_da * yc + pkm2 * c; + Scalar dqk_da = dqkm1_da * z - qkm1 - dqkm2_da * yc + qkm2 * c; + + if (qk != zero) { + Scalar ans_prev = ans; + ans = pk / qk; + + Scalar dans_da_prev = dans_da; + dans_da = (dpk_da - ans * dqk_da) / qk; + + if (mode == VALUE) { + if (numext::abs(ans_prev - ans) <= machep * numext::abs(ans)) { + break; + } + } else { + if (numext::abs(dans_da - dans_da_prev) <= machep) { + break; + } + } + } + + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + dpkm2_da = dpkm1_da; + dpkm1_da = dpk_da; + dqkm2_da = dqkm1_da; + dqkm1_da = dqk_da; + + if (numext::abs(pk) > big) { + pkm2 *= biginv; + pkm1 *= biginv; + qkm2 *= biginv; + qkm1 *= biginv; + + dpkm2_da *= biginv; + dpkm1_da *= biginv; + dqkm2_da *= biginv; + dqkm1_da *= biginv; + } + } + + /* Compute x**a * exp(-x) / gamma(a) */ + Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a); + Scalar dax_da = ax * dlogax_da; + + switch (mode) { + case VALUE: + return ans * ax; + case DERIVATIVE: + return ans * dax_da + dans_da * ax; + case SAMPLE_DERIVATIVE: + default: // this is needed to suppress clang warning + return -(dans_da + ans * dlogax_da) * x; + } + } +}; + +template <typename Scalar, IgammaComputationMode mode> +struct igamma_series_impl { + /* Computes igam(a, x) or its derivative (depending on the mode) + * using the series expansion of the incomplete Gamma function. + * + * Preconditions: + * x > 0 + * a > 0 + * !(x > 1 && x > a) + */ + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + const Scalar zero = 0; + const Scalar one = 1; + const Scalar machep = cephes_helper<Scalar>::machep(); + + Scalar ax = main_igamma_term<Scalar>(a, x); + + // This is independent of mode. If this value is zero, + // then the function value is zero. If the function value is zero, + // then we are in a neighborhood where the function value evalutes to zero, + // so the derivative is zero. + if (ax == zero) { + return zero; + } + + ax /= a; + + /* power series */ + Scalar r = a; + Scalar c = one; + Scalar ans = one; + + Scalar dc_da = zero; + Scalar dans_da = zero; + + for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) { + r += one; + Scalar term = x / r; + Scalar dterm_da = -x / (r * r); + dc_da = term * dc_da + dterm_da * c; + dans_da += dc_da; + c *= term; + ans += c; + + if (mode == VALUE) { + if (c <= machep * ans) { + break; + } + } else { + if (numext::abs(dc_da) <= machep * numext::abs(dans_da)) { + break; + } + } + } + + Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a + one); + Scalar dax_da = ax * dlogax_da; + + switch (mode) { + case VALUE: + return ans * ax; + case DERIVATIVE: + return ans * dax_da + dans_da * ax; + case SAMPLE_DERIVATIVE: + default: // this is needed to suppress clang warning + return -(dans_da + ans * dlogax_da) * x / a; + } + } +}; + #if !EIGEN_HAS_C99_MATH template <typename Scalar> @@ -466,8 +962,6 @@ struct igammac_impl { #else -template <typename Scalar> struct igamma_impl; // predeclare igamma_impl - template <typename Scalar> struct igammac_impl { EIGEN_DEVICE_FUNC @@ -535,93 +1029,15 @@ struct igammac_impl { return nan; } - if ((x < one) || (x < a)) { - /* The checks above ensure that we meet the preconditions for - * igamma_impl::Impl(), so call it, rather than igamma_impl::Run(). - * Calling Run() would also work, but in that case the compiler may not be - * able to prove that igammac_impl::Run and igamma_impl::Run are not - * mutually recursive. This leads to worse code, particularly on - * platforms like nvptx, where recursion is allowed only begrudgingly. - */ - return (one - igamma_impl<Scalar>::Impl(a, x)); - } - - return Impl(a, x); - } - - private: - /* igamma_impl calls igammac_impl::Impl. */ - friend struct igamma_impl<Scalar>; - - /* Actually computes igamc(a, x). - * - * Preconditions: - * a > 0 - * x >= 1 - * x >= a - */ - EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { - const Scalar zero = 0; - const Scalar one = 1; - const Scalar two = 2; - const Scalar machep = cephes_helper<Scalar>::machep(); - const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); - const Scalar big = cephes_helper<Scalar>::big(); - const Scalar biginv = cephes_helper<Scalar>::biginv(); - const Scalar inf = NumTraits<Scalar>::infinity(); - - Scalar ans, ax, c, yc, r, t, y, z; - Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; - - if (x == inf) return zero; // std::isinf crashes on CUDA - - /* Compute x**a * exp(-x) / gamma(a) */ - ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); - if (ax < -maxlog) { // underflow - return zero; + if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans + return nan; } - 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; - } + if ((x < one) || (x < a)) { + return (one - igamma_series_impl<Scalar, VALUE>::run(a, x)); } - return (ans * ax); + return igammac_cf_impl<Scalar, VALUE>::run(a, x); } }; @@ -631,15 +1047,10 @@ struct igammac_impl { * Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 * ************************************************************************************************/ -template <typename Scalar> -struct igamma_retval { - typedef Scalar type; -}; - #if !EIGEN_HAS_C99_MATH -template <typename Scalar> -struct igamma_impl { +template <typename Scalar, IgammaComputationMode mode> +struct igamma_generic_impl { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) { EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), @@ -650,69 +1061,17 @@ struct igamma_impl { #else -template <typename Scalar> -struct igamma_impl { +template <typename Scalar, IgammaComputationMode mode> +struct igamma_generic_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) + /* Depending on the mode, returns + * - VALUE: incomplete Gamma function igamma(a, x) + * - DERIVATIVE: derivative of incomplete Gamma function d/da igamma(a, x) + * - SAMPLE_DERIVATIVE: implicit derivative of a Gamma random variable + * x ~ Gamma(x | a, 1), dx/da = -1 / Gamma(x | a, 1) * d igamma(a, x) / dx * + * Derivatives are implemented by forward-mode differentiation. */ const Scalar zero = 0; const Scalar one = 1; @@ -724,67 +1083,167 @@ struct igamma_impl { return nan; } + if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans + return nan; + } + if ((x > one) && (x > a)) { - /* The checks above ensure that we meet the preconditions for - * igammac_impl::Impl(), so call it, rather than igammac_impl::Run(). - * Calling Run() would also work, but in that case the compiler may not be - * able to prove that igammac_impl::Run and igamma_impl::Run are not - * mutually recursive. This leads to worse code, particularly on - * platforms like nvptx, where recursion is allowed only begrudgingly. - */ - return (one - igammac_impl<Scalar>::Impl(a, x)); + Scalar ret = igammac_cf_impl<Scalar, mode>::run(a, x); + if (mode == VALUE) { + return one - ret; + } else { + return -ret; + } } - return Impl(a, x); + return igamma_series_impl<Scalar, mode>::run(a, x); } +}; + +#endif // EIGEN_HAS_C99_MATH - private: - /* igammac_impl calls igamma_impl::Impl. */ - friend struct igammac_impl<Scalar>; +template <typename Scalar> +struct igamma_retval { + typedef Scalar type; +}; - /* Actually computes igam(a, x). +template <typename Scalar> +struct igamma_impl : igamma_generic_impl<Scalar, VALUE> { + /* igam() + * Incomplete gamma integral. + * + * The CDF of Gamma(a, 1) random variable at the point x. + * + * Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample + * 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points. + * The ground truth is computed by mpmath. Mean absolute error: + * float: 1.26713e-05 + * double: 2.33606e-12 + * + * Cephes documentation below. + * + * 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 * - * Preconditions: - * x > 0 - * a > 0 - * !(x > 1 && x > a) */ - EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { - const Scalar zero = 0; - const Scalar one = 1; - const Scalar machep = cephes_helper<Scalar>::machep(); - const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); + /* + 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 + */ - Scalar ans, ax, c, r; + /* left tail of incomplete gamma function: + * + * inf. k + * a -x - x + * x e > ---------- + * - - + * k=0 | (a+k+1) + * + */ +}; - /* Compute x**a * exp(-x) / gamma(a) */ - ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); - if (ax < -maxlog) { - // underflow - return zero; - } - ax = numext::exp(ax); +template <typename Scalar> +struct igamma_der_a_retval : igamma_retval<Scalar> {}; - /* power series */ - r = a; - c = one; - ans = one; +template <typename Scalar> +struct igamma_der_a_impl : igamma_generic_impl<Scalar, DERIVATIVE> { + /* Derivative of the incomplete Gamma function with respect to a. + * + * Computes d/da igamma(a, x) by forward differentiation of the igamma code. + * + * Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample + * 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points. + * The ground truth is computed by mpmath. Mean absolute error: + * float: 6.17992e-07 + * double: 4.60453e-12 + * + * Reference: + * R. Moore. "Algorithm AS 187: Derivatives of the incomplete gamma + * integral". Journal of the Royal Statistical Society. 1982 + */ +}; - while (true) { - r += one; - c *= x/r; - ans += c; - if (c/ans <= machep) { - break; - } - } +template <typename Scalar> +struct gamma_sample_der_alpha_retval : igamma_retval<Scalar> {}; - return (ans * ax / a); - } +template <typename Scalar> +struct gamma_sample_der_alpha_impl + : igamma_generic_impl<Scalar, SAMPLE_DERIVATIVE> { + /* Derivative of a Gamma random variable sample with respect to alpha. + * + * Consider a sample of a Gamma random variable with the concentration + * parameter alpha: sample ~ Gamma(alpha, 1). The reparameterization + * derivative that we want to compute is dsample / dalpha = + * d igammainv(alpha, u) / dalpha, where u = igamma(alpha, sample). + * However, this formula is numerically unstable and expensive, so instead + * we use implicit differentiation: + * + * igamma(alpha, sample) = u, where u ~ Uniform(0, 1). + * Apply d / dalpha to both sides: + * d igamma(alpha, sample) / dalpha + * + d igamma(alpha, sample) / dsample * dsample/dalpha = 0 + * d igamma(alpha, sample) / dalpha + * + Gamma(sample | alpha, 1) dsample / dalpha = 0 + * dsample/dalpha = - (d igamma(alpha, sample) / dalpha) + * / Gamma(sample | alpha, 1) + * + * Here Gamma(sample | alpha, 1) is the PDF of the Gamma distribution + * (note that the derivative of the CDF w.r.t. sample is the PDF). + * See the reference below for more details. + * + * The derivative of igamma(alpha, sample) is computed by forward + * differentiation of the igamma code. Division by the Gamma PDF is performed + * in the same code, increasing the accuracy and speed due to cancellation + * of some terms. + * + * Accuracy estimation. For each alpha in [10^-2, 10^-1...10^3] we sample + * 50 Gamma random variables sample ~ Gamma(sample | alpha, 1), a total of 300 + * points. The ground truth is computed by mpmath. Mean absolute error: + * float: 2.1686e-06 + * double: 1.4774e-12 + * + * Reference: + * M. Figurnov, S. Mohamed, A. Mnih "Implicit Reparameterization Gradients". + * 2018 + */ }; -#endif // EIGEN_HAS_C99_MATH - /***************************************************************************** * Implementation of Riemann zeta function of two arguments, based on Cephes * *****************************************************************************/ @@ -944,7 +1403,12 @@ struct zeta_impl { { if(q == numext::floor(q)) { - return maxnum; + if (x == numext::floor(x) && long(x) % 2 == 0) { + return maxnum; + } + else { + return nan; + } } p = x; r = numext::floor(p); @@ -1020,11 +1484,11 @@ struct polygamma_impl { Scalar nplus = n + one; const Scalar nan = NumTraits<Scalar>::quiet_NaN(); - // Check that n is an integer - if (numext::floor(n) != n) { + // Check that n is a non-negative integer + if (numext::floor(n) != n || n < zero) { return nan; } - // Just return the digamma function for n = 1 + // Just return the digamma function for n = 0 else if (n == zero) { return digamma_impl<Scalar>::run(x); } @@ -1392,7 +1856,7 @@ struct betainc_helper<double> { if ((a + b) < maxgam && numext::abs(u) < maxlog) { t = gamma(a + b) / (gamma(a) * gamma(b)); s = s * t * pow(x, a); - } else { + } */ t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) - lgamma_impl<double>::run(b) + u + numext::log(s); @@ -1540,12 +2004,30 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) } template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(ndtri, Scalar) + ndtri(const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(ndtri, Scalar)::run(x); +} + +template <typename Scalar> EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar) igamma(const Scalar& a, const Scalar& x) { return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x); } template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma_der_a, Scalar) + igamma_der_a(const Scalar& a, const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(igamma_der_a, Scalar)::run(a, x); +} + +template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(gamma_sample_der_alpha, Scalar) + gamma_sample_der_alpha(const Scalar& a, const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(gamma_sample_der_alpha, Scalar)::run(a, x); +} + +template <typename Scalar> EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar) igammac(const Scalar& a, const Scalar& x) { return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x); @@ -1558,8 +2040,6 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(betainc, Scalar) } } // end namespace numext - - } // end namespace Eigen #endif // EIGEN_SPECIAL_FUNCTIONS_H |