GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-10 18:59:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix edge-case in zeta for large inputs.

Browse Source
This commit is contained in:
Antonio Sánchez 2022-03-08 21:21:20 +00:00 committed by Rasmus Munk Larsen
parent cd2ba9d03e
commit 9296bb4b93
2 changed files with 13 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1388,7 +1388,7 @@ struct zeta_impl {
}; };
const Scalar maxnum = NumTraits<Scalar>::infinity(); const Scalar maxnum = NumTraits<Scalar>::infinity();
const Scalar zero = 0.0, half = 0.5, one = 1.0; const Scalar zero = Scalar(0.0), half = Scalar(0.5), one = Scalar(1.0);
const Scalar machep = cephes_helper<Scalar>::machep(); const Scalar machep = cephes_helper<Scalar>::machep();
const Scalar nan = NumTraits<Scalar>::quiet_NaN(); const Scalar nan = NumTraits<Scalar>::quiet_NaN();
@ -1430,11 +1430,19 @@ struct zeta_impl {
return s; return s;
} }
// If b is zero, then the tail sum will also end up being zero.
// Exiting early here can prevent NaNs for some large inputs, where
// the tail sum computed below has term `a` which can overflow to `inf`.
if (numext::equal_strict(b, zero)) {
return s;
}
w = a; w = a;
s += b*w/(x-one); s += b*w/(x-one);
s -= half * b; s -= half * b;
a = one; a = one;
k = zero; k = zero;
for( i=0; i<12; i++ ) for( i=0; i<12; i++ )
{ {
a *= x + k; a *= x + k;

8
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(10), q(10), res(10), ref(10); ArrayType x(11), q(11), res(11), ref(11);
x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9, 2, 3, 4; x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9, 2, 3, 4, 2000;
q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345, -1, -2, -3; q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345, -1, -2, -3, 2000;
ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan, plusinf, nan, plusinf; ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan, plusinf, nan, plusinf, 0;
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); );