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

Fixed suggestions by Eugene Brevdo.

Browse Source
This commit is contained in:
Till Hoffmann 2016-04-01 17:51:39 +01:00
parent 57239f4a81
commit 3cb0a237c1
2 changed files with 52 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
Eigen/src/Core/SpecialFunctions.h
View File

@ -806,10 +806,9 @@ struct zeta_impl {
*/
int i;
/*double a, b, k, s, t, w;*/
Scalar p, r, a, b, k, s, t, w;
const double A[] = {
const Scalar A[] = {
12.0,
-720.0,
30240.0,
@ -829,12 +828,10 @@ struct zeta_impl {
const Scalar machep = igamma_helper<Scalar>::machep();
if( x == one )
return maxnum; //goto retinf;
return maxnum;
if( x < one )
{
// domerr:
// mtherr( "zeta", DOMAIN );
return zero;
}
@ -842,64 +839,53 @@ struct zeta_impl {
{
if(q == numext::floor(q))
{
// mtherr( "zeta", SING );
// retinf:
return maxnum;
}
p = x;
r = numext::floor(p);
// if( x != floor(x) )
// goto domerr; /* because q^-x not defined */
if (p != r)
return zero;
}
/* Euler-Maclaurin summation formula */
/*
if( x < 25.0 )
/* Permit negative q but continue sum until n+q > +9 .
* This case should be handled by a reflection formula.
* If q<0 and x is an integer, there is a relation to
* the polygamma function.
*/
s = numext::pow( q, -x );
a = q;
i = 0;
b = zero;
while( (i < 9) || (a <= Scalar(9)) )
{
/* Permit negative q but continue sum until n+q > +9 .
* This case should be handled by a reflection formula.
* If q<0 and x is an integer, there is a relation to
* the polygamma function.
*/
s = numext::pow( q, -x );
a = q;
i = 0;
b = zero;
while( (i < 9) || (a <= Scalar(9.0)) )
{
i += 1;
a += one;
b = numext::pow( a, -x );
s += b;
if( numext::abs(b/s) < machep )
return s; // goto done;
}
w = a;
s += b*w/(x-one);
s -= half * b;
a = one;
k = zero;
for( i=0; i<12; i++ )
{
a *= x + k;
b /= w;
t = a*b/A[i];
s = s + t;
t = numext::abs(t/s);
if( t < machep )
return s; // goto done;
k += one;
a *= x + k;
b /= w;
k += one;
}
// done:
return(s);
}
i += 1;
a += one;
b = numext::pow( a, -x );
s += b;
if( numext::abs(b/s) < machep )
return s;
}
w = a;
s += b*w/(x-one);
s -= half * b;
a = one;
k = zero;
for( i=0; i<12; i++ )
{
a *= x + k;
b /= w;
t = a*b/A[i];
s = s + t;
t = numext::abs(t/s);
if( t < machep )
return s;
k += one;
a *= x + k;
b /= w;
k += one;
}
return s;
}
};

14
test/array.cpp
View File

@ -332,10 +332,22 @@ template<typename ArrayType> void array_real(const ArrayType& m)
VERIFY_IS_EQUAL(numext::zeta(Scalar(1), Scalar(1.2345)), // The second scalar does not matter
std::numeric_limits<RealScalar>::infinity());
// Check the polygamma against scipy.special.polygamma
// Check the polygamma against scipy.special.polygamma examples
VERIFY_IS_APPROX(numext::polygamma(Scalar(1), Scalar(2)), RealScalar(0.644934066848));
VERIFY_IS_APPROX(numext::polygamma(Scalar(1), Scalar(3)), RealScalar(0.394934066848));
VERIFY_IS_APPROX(numext::polygamma(Scalar(1), Scalar(25.5)), RealScalar(0.0399946696496));
// Check the polygamma function over a larger range of values
VERIFY_IS_APPROX(numext::polygamma(Scalar(17), Scalar(4.7)), RealScalar(293.334565435));
VERIFY_IS_APPROX(numext::polygamma(Scalar(31), Scalar(11.8)), RealScalar(0.445487887616));
VERIFY_IS_APPROX(numext::polygamma(Scalar(28), Scalar(17.7)), RealScalar(-2.47810300902e-07));
VERIFY_IS_APPROX(numext::polygamma(Scalar(8), Scalar(30.2)), RealScalar(-8.29668781082e-09));
/* The following tests only pass for doubles because floats cannot handle the large values of
the gamma function.
VERIFY_IS_APPROX(numext::polygamma(Scalar(42), Scalar(15.8)), RealScalar(-0.434562276666));
VERIFY_IS_APPROX(numext::polygamma(Scalar(147), Scalar(54.1)), RealScalar(0.567742190178));
VERIFY_IS_APPROX(numext::polygamma(Scalar(170), Scalar(64)), RealScalar(-0.0108615497927));
*/
{
// Test various propreties of igamma & igammac. These are normalized