diff options
| -rw-r--r-- | Eigen/src/Core/SpecialFunctions.h | 119 | ||||
| -rw-r--r-- | test/array.cpp | 64 |
2 files changed, 119 insertions, 64 deletions
diff --git a/Eigen/src/Core/SpecialFunctions.h b/Eigen/src/Core/SpecialFunctions.h index ff2146afc..4a61325d4 100644 --- a/Eigen/src/Core/SpecialFunctions.h +++ b/Eigen/src/Core/SpecialFunctions.h @@ -296,7 +296,8 @@ struct digamma_impl { if (x <= zero) { negative = one; q = x; - p = ::floor(q); + using std::floor; + p = floor(q); if (p == q) { return maxnum; } @@ -309,7 +310,8 @@ struct digamma_impl { p += one; nz = q - p; } - nz = m_pi / ::tan(m_pi * nz); + using std::tan; + nz = m_pi / tan(m_pi * nz); } else { nz = zero; @@ -327,7 +329,8 @@ struct digamma_impl { y = digamma_impl_maybe_poly<Scalar>::run(s); - y = ::log(s) - (half / s) - y - w; + using std::log; + y = log(s) - (half / s) - y - w; return (negative) ? y - nz : y; } @@ -427,6 +430,39 @@ struct igammac_impl { template <typename Scalar> struct igamma_impl; // predeclare igamma_impl template <typename Scalar> +struct igamma_helper { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + static Scalar machep() { assert(false && "machep not supported for this type"); return 0.0; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + static Scalar big() { assert(false && "big not supported for this type"); return 0.0; } +}; + +template <> +struct igamma_helper<float> { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + static float machep() { + return NumTraits<float>::epsilon() / 2; // 1.0 - machep == 1.0 + } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + static float big() { + // use epsneg (1.0 - epsneg == 1.0) + return 1.0 / (NumTraits<float>::epsilon() / 2); + } +}; + +template <> +struct igamma_helper<double> { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + static double machep() { + return NumTraits<double>::epsilon() / 2; // 1.0 - machep == 1.0 + } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + static double big() { + return 1.0 / NumTraits<double>::epsilon(); + } +}; + +template <typename Scalar> struct igammac_impl { EIGEN_DEVICE_FUNC static Scalar run(Scalar a, Scalar x) { @@ -487,26 +523,35 @@ struct igammac_impl { const Scalar zero = 0; const Scalar one = 1; const Scalar two = 2; - const Scalar machep = NumTraits<Scalar>::epsilon(); + const Scalar machep = igamma_helper<Scalar>::machep(); const Scalar maxlog = ::log(NumTraits<Scalar>::highest()); - const Scalar big = one / machep; + const Scalar big = igamma_helper<Scalar>::big(); + const Scalar biginv = 1 / big; + const Scalar nan = NumTraits<Scalar>::quiet_NaN(); Scalar ans, ax, c, yc, r, t, y, z; Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; - if ((x <= zero) || ( a <= zero)) { - return one; + if ((x < zero) || ( a <= zero)) { + // domain error + return nan; } if ((x < one) || (x < a)) { return (one - igamma_impl<Scalar>::run(a, x)); } - ax = a * ::log(x) - x - lgamma_impl<Scalar>::run(a); - if( ax < -maxlog ) { // underflow + using std::isinf; + if ((isinf)(x)) return zero; + + /* Compute x**a * exp(-x) / gamma(a) */ + using std::log; + ax = a * log(x) - x - lgamma_impl<Scalar>::run(a); + if (ax < -maxlog) { // underflow return zero; } - ax = ::exp(ax); + using std::exp; + ax = exp(ax); // continued fraction y = one - a; @@ -516,35 +561,36 @@ struct igammac_impl { qkm2 = x; pkm1 = x + one; qkm1 = z * x; - ans = pkm1/qkm1; + ans = pkm1 / qkm1; + using std::abs; do { 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 = ::abs( (ans - r)/r ); + pk = pkm1 * z - pkm2 * yc; + qk = qkm1 * z - qkm2 * yc; + if (qk != zero) { + r = pk / qk; + t = abs((ans - r) / r); ans = r; - } else { + } else { t = one; } pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; - if (::abs(pk) > big) { - pkm2 *= machep; - pkm1 *= machep; - qkm2 *= machep; - qkm1 *= machep; + if (abs(pk) > big) { + pkm2 *= biginv; + pkm1 *= biginv; + qkm2 *= biginv; + qkm1 *= biginv; } - } while( t > machep ); + } while (t > machep); - return ( ans * ax ); + return (ans * ax); } }; @@ -639,26 +685,31 @@ struct igamma_impl { */ const Scalar zero = 0; const Scalar one = 1; - const Scalar machep = NumTraits<Scalar>::epsilon(); + const Scalar machep = igamma_helper<Scalar>::machep(); const Scalar maxlog = ::log(NumTraits<Scalar>::highest()); + const Scalar nan = NumTraits<Scalar>::quiet_NaN(); double ans, ax, c, r; - if( (x <= zero) || ( a <= zero) ) { - return zero; + if (x == zero) return zero; + + if ((x < zero) || ( a <= zero)) { // domain error + return nan; } - if( (x > one) && (x > a ) ) { - return (one - igammac_impl<Scalar>::run(a,x)); + if ((x > one) && (x > a)) { + return (one - igammac_impl<Scalar>::run(a, x)); } /* Compute x**a * exp(-x) / gamma(a) */ - ax = a * ::log(x) - x - lgamma_impl<Scalar>::run(a); - if( ax < -maxlog ) { + using std::log; + ax = a * log(x) - x - lgamma_impl<Scalar>::run(a); + if (ax < -maxlog) { // underflow return zero; } - ax = ::exp(ax); + using std::exp; + ax = exp(ax); /* power series */ r = a; @@ -669,9 +720,9 @@ struct igamma_impl { r += one; c *= x/r; ans += c; - } while( c/ans > machep ); + } while (c/ans > machep); - return( ans * ax/a ); + return (ans * ax / a); } }; diff --git a/test/array.cpp b/test/array.cpp index a37874cc2..c61bfc8ed 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -331,41 +331,45 @@ template<typename ArrayType> void array_real(const ArrayType& m) VERIFY_IS_EQUAL(numext::digamma(Scalar(-1)), std::numeric_limits<RealScalar>::infinity()); - Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(10000.5)}; - Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(10000.5)}; + 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<Scalar>::quiet_NaN(); - Scalar igamma_s[][6] = { - {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, - {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.5013297751014064}}; - Scalar igammac_s[][6] = { - {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, - {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.49867022490946517}}; + 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) { - //std::cout << numext::igamma(a_s[i], x_s[j]) << " vs. " << igamma_s[i][j] << std::endl; - //std::cout << numext::igammac(a_s[i], x_s[j]) << " c.vs. " << - //igammac_s[i][j] << std::endl; - std::cout << a_s[i] << ", " << x_s[j] << std::endl; - VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]); - VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][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]); + } } } } |