boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

[polygamma] Fix spurious overflow in transition zone, update cot derivative to work to arbitrary level.

Browse Source
This commit is contained in:
jzmaddock
2014-11-03 11:50:10 +00:00
parent e70f53e670
commit eccec791c7
3 changed files with 112 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

109
include/boost/math/special_functions/detail/polygamma.hpp
View File

@@ -134,14 +134,25 @@
T sum0(0);
T z_plus_k_pow_minus_m_minus_one(0);
// Forward recursion to larger x:
for(int k = 1; k <= iter; ++k)
// Forward recursion to larger x, need to check for overflow first though:
if(log(z + iter) * minus_m_minus_one > -tools::log_max_value<T>())
{
z_plus_k_pow_minus_m_minus_one = pow(z, minus_m_minus_one);
sum0 += z_plus_k_pow_minus_m_minus_one;
z += 1;
for(int k = 1; k <= iter; ++k)
{
z_plus_k_pow_minus_m_minus_one = pow(z, minus_m_minus_one);
sum0 += z_plus_k_pow_minus_m_minus_one;
z += 1;
}
sum0 *= boost::math::factorial<T>(n);
}
{
for(int k = 1; k <= iter; ++k)
{
T log_term = log(z) * minus_m_minus_one + boost::math::lgamma<T>(n + 1, pol);
sum0 += exp(log_term);
z += 1;
}
}
sum0 *= boost::math::factorial<T>(n);
if((n - 1) & 1)
sum0 = -sum0;
@@ -232,6 +243,18 @@
return n & 1 ? sum : -sum;
}
//
// Helper function which figures out which slot our coefficient is in
// given an angle multiplier for the cosine term of power:
//
template <class Table>
typename Table::value_type::reference dereference_table(Table& table, unsigned row, unsigned power)
{
return table[row][power / 2];
}
template <class T, class Policy>
T poly_cot_pi(int n, T x, T xc, const Policy& pol, const char* function)
{
@@ -252,7 +275,7 @@
// of storage space, this method has no better accuracy than recursion
// on x to x > 0 when computing polygamma :-(
//
T s = x < xc ? boost::math::sin_pi(x) : boost::math::sin_pi(xc);
T s = fabs(x) < fabs(xc) ? boost::math::sin_pi(x) : boost::math::sin_pi(xc);
switch(n)
{
case 1:
@@ -447,11 +470,67 @@
+ 2097130 * boost::math::cos_pi(18 * x)
+ boost::math::cos_pi(20 * x)) / boost::math::pow<22>(s);
}
/*
*/
}
return policies::raise_domain_error<T>(function, "Derivative of cotangent not implemented at n = %1%", n, pol);
//
// We'll have to compute the corefficients up to n:
//
static std::vector<std::vector<T> > table(1, std::vector<T>(1, T(-1)));
int index = n - 1;
if(index >= (int)table.size())
{
//
// We need to compute new coefficients for the cosine terms to the derivative.
// The following code follows immediately from differentiating
//
// C * cos(power * x) * csc^n(power * x)
//
for(int i = table.size(); i <= index; ++i)
{
table.push_back(std::vector<T>((i + 2) / 2, T(0)));
for(int power = (i & 1 ? 0 : 1); power < i; power += 2)
{
dereference_table(table, i, std::abs(1 - power)) += -(power + i + 1) * dereference_table(table, i - 1, power) / 2;
dereference_table(table, i, power + 1) += -(i + 1 - power) * dereference_table(table, i - 1, power) / 2;
}
//
// The coefficients as calculated above grow so large so fast that we scale them all
// by n! And since current order = i+1 just divide each row by that as we create it:
//
for(unsigned j = 0; j < table[i].size(); ++j)
table[i][j] /= (i + 1);
}
}
T sum = 0;
int power = n & 1 ? 0 : 1;
//
// Compute the sum of the cosine terms:
//
for(unsigned j = 0; j < table[index].size(); ++j)
{
sum += table[index][j] * boost::math::cos_pi(x * power);
power += 2;
}
if(sum == 0)
return sum;
//
// The remaining terms are computed using logs since the powers and factorials
// get real large real quick:
//
T power_terms = n * log(boost::math::constants::pi<T>());
if(s == 0)
return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol);
power_terms -= log(fabs(s)) * (n + 1);
power_terms += boost::math::lgamma(T(n + 1));
power_terms += log(fabs(sum));
if(power_terms > boost::math::tools::log_max_value<T>())
return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol);
return exp(power_terms) * ((s < 0) && ((n + 1) & 1) ? -1 : 1) * boost::math::sign(sum);
}
template<class T, class Policy>
@@ -475,15 +554,16 @@
else
return policies::raise_pole_error<T>(function, "Evaluation at negative integer %1%", x, pol);
}
if(n < 22)
{
//if(n < 22)
//{
//
// We have tabulated the derivatives of cot(x) up to the 9th derivative, which
// allows us to use: http://functions.wolfram.com/06.15.16.0001.01
T z = 1 - x;
T result = polygamma_imp(n, z, pol) + constants::pi<T, Policy>() * poly_cot_pi(n, z, x, pol, function);
return n & 1 ? T(-result) : result;
}
//}
#if 0
//
// Try http://functions.wolfram.com/06.15.16.0007.01
//
@@ -500,6 +580,7 @@
if(n & 1)
sum = -sum;
return polygamma_imp(n, z, pol) - sum;
#endif
}
//
// Limit for use of small-x-series is chosen

7
test/test_polygamma.cpp
View File

@@ -41,6 +41,13 @@ void expected_results()
largest_type, // test type(s)
".*negative.*", // test data group
".*", 1400, 500); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
largest_type, // test type(s)
".*bug cases.*", // test data group
".*", 100000, 50000); // test function
add_expected_result(
".*", // compiler
".*", // stdlib

12
test/test_polygamma.hpp
View File

@@ -160,8 +160,6 @@ void test_polygamma(T, const char* name)
{ SC_(2.0), SC_(-9.5), SC_(-0.0099751442477151692853059194570941025) }, { SC_(4.0), SC_(-9.5), SC_(-0.00059506011900940675655749713967447346) }, { SC_(6.0), SC_(-9.5), SC_(-0.00011794286977626608581527674104044053) }, { SC_(8.0), SC_(-9.5), SC_(-0.000048934615584055214532361558113243801) }, { SC_(10.0), SC_(-9.5), SC_(-0.000034696555222805555969152083201249795) }, { SC_(12.0), SC_(-9.5), SC_(-0.000037470635416758472254487967117333555) }, { SC_(14.0), SC_(-9.5), SC_(-0.000057218576281198884425118075685027774) }, { SC_(16.0), SC_(-9.5), SC_(-0.00011728023376485851827598955099805232) }, { SC_(18.0), SC_(-9.5), SC_(-0.00031049110045758006635527458576736345) }, { SC_(20.0), SC_(-9.5), SC_(-0.0010307637762451416598018081796345932) },
{ SC_(2.0), SC_(-9.5367431640625e-7), SC_(2.3058430092136939495958800005662742e18) }, { SC_(2.0), SC_(-5.960464477539063e-8), SC_(9.4447329657392904273895958858066118e21) }, { SC_(2.0), SC_(-3.725290298461914e-9), SC_(3.8685626227668133590597629595886169e25) }, { SC_(2.0), SC_(-2.3283064365386963e-10), SC_(1.5845632502852867518708790066959589e29) }, { SC_(2.0), SC_(-1.4551915228366852e-11), SC_(6.4903710731685345356631204115250960e32) }, { SC_(2.0), SC_(-9.094947017729282e-13), SC_(2.6584559915698317458076141205606891e36) }, { SC_(2.0), SC_(-5.684341886080802e-14), SC_(1.0889035741470030830827987437816583e40) }, { SC_(2.0), SC_(-3.552713678800501e-15), SC_(4.4601490397061246283071436545296723e43) }, { SC_(2.0), SC_(-2.220446049250313e-16), SC_(1.8268770466636286477546060408953538e47) }, { SC_(2.0), SC_(-1.3877787807814457e-17), SC_(7.4828883831342229412028663435073691e50) }, { SC_(2.0), SC_(-8.673617379884035e-19), SC_(3.0649910817317777167166940543006184e54) }, { SC_(2.0), SC_(-5.421010862427522e-20), SC_(1.2554203470773361527671578846415333e58) }, { SC_(2.0), SC_(-3.3881317890172014e-21), SC_(5.1422017416287688817342786954917203e61) },
{ SC_(3.0), SC_(-9.5367431640625e-7), SC_(7.2535549176877750482370624939631357e24) }, { SC_(3.0), SC_(-5.960464477539063e-8), SC_(4.7536897508558602556126370202249394e29) }, { SC_(3.0), SC_(-3.725290298461914e-9), SC_(3.1153781151208965771182977975320582e34) }, { SC_(3.0), SC_(-2.3283064365386963e-10), SC_(2.0416942015256307807802476445906093e39) }, { SC_(3.0), SC_(-1.4551915228366852e-11), SC_(1.3380447119118373884921430963589017e44) }, { SC_(3.0), SC_(-9.094947017729282e-13), SC_(8.7690098239854175092221089962976981e48) }, { SC_(3.0), SC_(-5.684341886080802e-14), SC_(5.7468582782470832188438013518136594e53) }, { SC_(3.0), SC_(-3.552713678800501e-15), SC_(3.7662610412320084583014736539245998e58) }, { SC_(3.0), SC_(-2.220446049250313e-16), SC_(2.4682568359818090632324537738360258e63) }, { SC_(3.0), SC_(-1.3877787807814457e-17), SC_(1.6175968000290383876800209052211778e68) }, { SC_(3.0), SC_(-8.673617379884035e-19), SC_(1.0601082388670305977499785004457511e73) }, { SC_(3.0), SC_(-5.421010862427522e-20), SC_(6.9475253542389717254142591005212745e77) }, { SC_(3.0), SC_(-3.3881317890172014e-21), SC_(4.5531302161540525099674888441176224e82) },
//{ SC_(20.0), SC_(-9.5), SC_(-0.0010307637762451416598018081796345932)}, { SC_(21.0), SC_(-9.5), SC_(4.2858188624279670465725116159418790e26) }, { SC_(22.0), SC_(-9.5), SC_(-0.0041914420015886426448664611125724338) }, { SC_(23.0), SC_(-9.5), SC_(8.6744973773084910367753904944787155e29) }, { SC_(24.0), SC_(-9.5), SC_(-0.020482527998674199369359987617445420) }, { SC_(25.0), SC_(-9.5), SC_(2.0818793705474855280713094147422278e33) }, { SC_(26.0), SC_(-9.5), SC_(-0.11840299831136879082790399023048278) }, { SC_(27.0), SC_(-9.5), SC_(5.8459172724952950454469355072783445e36) }, { SC_(28.0), SC_(-9.5), SC_(-0.79896867888629450211348696225656872) }, { SC_(29.0), SC_(-9.5), SC_(1.8987539301063980537054871851063348e40) }, { SC_(30.0), SC_(-9.5), SC_(-6.2224465669723272001827279817965087) }
} };
do_test_polygamma<T>(neg_data, name, "Mathematica Data - negative arguments");
@@ -179,5 +177,15 @@ void test_polygamma(T, const char* name)
{ SC_(35.0), SC_(0.1250000000), SC_(3.3532982327901451846973629635910627e72) }, { SC_(35.0), SC_(0.06250000000), SC_(2.3043689989709229438987285737704404e83) }, { SC_(35.0), SC_(0.03125000000), SC_(1.5835503181594194718369731136519624e94) }, { SC_(35.0), SC_(0.01562500000), SC_(1.0882074924904162473416628106826351e105) }, { SC_(35.0), SC_(0.007812500000), SC_(7.4781049464136054012151819362261716e115) }, { SC_(35.0), SC_(0.003906250000), SC_(5.1389145889443628300269064119650184e126) }, { SC_(35.0), SC_(0.001953125000), SC_(3.5314352154325314429637711743107931e137) }, { SC_(35.0), SC_(0.0009765625000), SC_(2.4267838013160699267233738410387795e148) }
} };
do_test_polygamma<T>(small_data, name, "Mathematica Data - small arguments");
boost::array<boost::array<value_type, 3>, 18> bug_cases =
{ {
{ SC_(171.0), SC_(2.0), SC_(2.073093314165313149880140394410e257) }, { SC_(171.0), SC_(5.0), SC_(7.42911976071332889749264626321716781e188) },
{ SC_(166.0), SC_(2.0), SC_(-4.8129498903508823293044351695484095e247) }, { SC_(166.0), SC_(3.0), SC_(-1.8843912448604502196243093626013895e218) },
{ SC_(171.0), SC_(23.0), SC_(7.53143916217078889612817829861181739e74) }, { SC_(168.0), SC_(150.0), SC_(-6.5266062780306068333215312257902920e-66) },
{ SC_(169.0), SC_(202.0), SC_(9.2734049986021958613510169328055599e-88) },
{ SC_(20.0), SC_(-9.5), SC_(-0.0010307637762451416598018081796345932)}, { SC_(21.0), SC_(-9.5), SC_(4.2858188624279670465725116159418790e26) }, { SC_(22.0), SC_(-9.5), SC_(-0.0041914420015886426448664611125724338) }, { SC_(23.0), SC_(-9.5), SC_(8.6744973773084910367753904944787155e29) }, { SC_(24.0), SC_(-9.5), SC_(-0.020482527998674199369359987617445420) }, { SC_(25.0), SC_(-9.5), SC_(2.0818793705474855280713094147422278e33) }, { SC_(26.0), SC_(-9.5), SC_(-0.11840299831136879082790399023048278) }, { SC_(27.0), SC_(-9.5), SC_(5.8459172724952950454469355072783445e36) }, { SC_(28.0), SC_(-9.5), SC_(-0.79896867888629450211348696225656872) }, { SC_(29.0), SC_(-9.5), SC_(1.8987539301063980537054871851063348e40) }, { SC_(30.0), SC_(-9.5), SC_(-6.2224465669723272001827279817965087) }
} };
do_test_polygamma<T>(bug_cases, name, "Mathematica Data - Large orders and other bug cases");
}