mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

Return nan at poles of polygamma, digamma, and zeta if limit is not defined

Browse Source
This commit is contained in:
frgossen 2021-02-19 16:35:11 +00:00 committed by Antonio Sánchez
parent 7f09d3487d
commit 33e0af0130
2 changed files with 22 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
View File

@ -241,7 +241,7 @@ struct digamma_impl {
Scalar p, q, nz, s, w, y; Scalar p, q, nz, s, w, y;
bool negative = false; 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 m_pi = Scalar(EIGEN_PI);
const Scalar zero = Scalar(0); const Scalar zero = Scalar(0);
@ -254,7 +254,7 @@ struct digamma_impl {
q = x; q = x;
p = numext::floor(q); p = numext::floor(q);
if (p == q) { if (p == q) {
return maxnum; return nan;
} }
/* Remove the zeros of tan(m_pi x) /* Remove the zeros of tan(m_pi x)
* by subtracting the nearest integer from x * by subtracting the nearest integer from x
@ -1403,7 +1403,12 @@ struct zeta_impl {
{ {
if(q == numext::floor(q)) if(q == numext::floor(q))
{ {
return maxnum; if (x == numext::floor(x) && long(x) % 2 == 0) {
return maxnum;
}
else {
return nan;
}
} }
p = x; p = x;
r = numext::floor(p); r = numext::floor(p);
@ -1479,11 +1484,11 @@ struct polygamma_impl {
Scalar nplus = n + one; Scalar nplus = n + one;
const Scalar nan = NumTraits<Scalar>::quiet_NaN(); const Scalar nan = NumTraits<Scalar>::quiet_NaN();
// Check that n is an integer // Check that n is a non-negative integer
if (numext::floor(n) != n) { if (numext::floor(n) != n || n < zero) {
return nan; return nan;
} }
// Just return the digamma function for n = 1 // Just return the digamma function for n = 0
else if (n == zero) { else if (n == zero) {
return digamma_impl<Scalar>::run(x); return digamma_impl<Scalar>::run(x);
} }

22
unsupported/test/special_functions.cpp
View File

@ -191,10 +191,10 @@ template<typename ArrayType> void array_special_functions()
// Check the zeta function against scipy.special.zeta // Check the zeta function against scipy.special.zeta
{ {
ArrayType x(7), q(7), res(7), ref(7); ArrayType x(10), q(10), res(10), ref(10);
x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9; x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9, 2, 3, 4;
q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345; q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345, -1, -2, -3;
ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan; ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan, plusinf, nan, plusinf;
CALL_SUBTEST( verify_component_wise(ref, ref); ); CALL_SUBTEST( verify_component_wise(ref, ref); );
CALL_SUBTEST( res = x.zeta(q); verify_component_wise(res, ref); ); CALL_SUBTEST( res = x.zeta(q); verify_component_wise(res, ref); );
CALL_SUBTEST( res = zeta(x,q); verify_component_wise(res, ref); ); CALL_SUBTEST( res = zeta(x,q); verify_component_wise(res, ref); );
@ -202,9 +202,9 @@ template<typename ArrayType> void array_special_functions()
// digamma // digamma
{ {
ArrayType x(7), res(7), ref(7); ArrayType x(9), res(9), ref(9);
x << 1, 1.5, 4, -10.5, 10000.5, 0, -1; x << 1, 1.5, 4, -10.5, 10000.5, 0, -1, -2, -3;
ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, plusinf, plusinf; ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, nan, nan, nan, nan;
CALL_SUBTEST( verify_component_wise(ref, ref); ); CALL_SUBTEST( verify_component_wise(ref, ref); );
CALL_SUBTEST( res = x.digamma(); verify_component_wise(res, ref); ); CALL_SUBTEST( res = x.digamma(); verify_component_wise(res, ref); );
@ -213,10 +213,10 @@ template<typename ArrayType> void array_special_functions()
#if EIGEN_HAS_C99_MATH #if EIGEN_HAS_C99_MATH
{ {
ArrayType n(11), x(11), res(11), ref(11); ArrayType n(16), x(16), res(16), ref(16);
n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170; n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170, -1, 0, 1, 2, 3;
x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64; x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64, -1, -2, -3, -4, -5;
ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927; ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927, nan, nan, plusinf, nan, plusinf;
CALL_SUBTEST( verify_component_wise(ref, ref); ); CALL_SUBTEST( verify_component_wise(ref, ref); );
if(sizeof(RealScalar)>=8) { // double if(sizeof(RealScalar)>=8) { // double