mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-09-12 09:23:12 +08:00
Make igamma and igammac work correctly.
This required replacing ::abs with std::abs. Modified some unit tests.
This commit is contained in:
Eugene Brevdo
parent
7ea35bfa1c
commit
0b9e0abc96
@ -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;
|
||||
}
|
||||
@ -426,6 +429,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
|
||||
@ -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);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@ -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]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user