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 issues with small negative arguments for x.

Browse Source 
Tidy up a couple of cosmetic coding issues.
This commit is contained in:
jzmaddock
2014-10-31 18:51:14 +00:00
parent 8438b8a84b
commit e70f53e670
2 changed files with 11 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -69,7 +69,7 @@
// set the initial values of part_term, term and sum via logs:
part_term = boost::math::lgamma(n, pol) - (n + 1) * log(x);
sum = exp(part_term + log(n + 2 * x) - boost::math::constants::ln_two<T>());
part_term += log(n * (n + 1)) - boost::math::constants::ln_two<T>() - log(x);
part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two<T>() - log(x);
part_term = exp(part_term);
}
else
@@ -123,9 +123,8 @@
// Use N = (0.4 * digits) + (4 * n) for target value for x:
BOOST_MATH_STD_USING
const int d4d = static_cast<boost::int32_t>(0.4F * policies::digits_base10<T, Policy>());
const int N4dn = static_cast<boost::int32_t>(d4d + (4 * n));
const int N = static_cast<boost::int32_t>((std::min)(N4dn, (std::numeric_limits<int>::max)()));
const int d4d = static_cast<int>(0.4F * policies::digits_base10<T, Policy>());
const int N = d4d + (4 * n);
const int m = n;
const int iter = N - itrunc(x);
@@ -234,7 +233,7 @@
}
template <class T, class Policy>
T poly_cot_pi(int n, T x, const Policy& pol, const char* function)
T poly_cot_pi(int n, T x, T xc, const Policy& pol, const char* function)
{
// Return n'th derivative of cot(pi*x) at x, these are simply
// tabulated for up to n = 9, beyond that it is possible to
@@ -253,7 +252,7 @@
// of storage space, this method has no better accuracy than recursion
// on x to x > 0 when computing polygamma :-(
//
T s = boost::math::sin_pi(x);
T s = x < xc ? boost::math::sin_pi(x) : boost::math::sin_pi(xc);
switch(n)
{
case 1:
@@ -482,7 +481,7 @@
// 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, pol, function);
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;
}
//

7
test/test_polygamma.hpp
View File

@@ -137,7 +137,7 @@ void test_polygamma(T, const char* name)
} };
do_test_polygamma<T>(big_data, name, "Mathematica Data - large arguments");
boost::array<boost::array<value_type, 3>, 373> neg_data =
boost::array<boost::array<value_type, 3>, 399> neg_data =
{ {
{ 1, SC_(-12.750), SC_(19.663772856722737612034697464751605) }, { 1, SC_(-12.250), SC_(19.660817549236368273654684043826967) }, { 1, SC_(-11.750), SC_(19.657621376522814505537196503582823) }, { 1, SC_(-11.250), SC_(19.654153659190554029589711115880278) }, { 1, SC_(-10.750), SC_(19.650378280099093364749509767503149) }, { 1, SC_(-10.250), SC_(19.646252424622652795021809881312377) }, { 1, SC_(-9.7500), SC_(19.641724953976865133273035997897957) }, { 1, SC_(-9.2500), SC_(19.636734280660725370869519577921538) }, { 1, SC_(-8.7500), SC_(19.631205558842085383108670448917024) }, { 1, SC_(-8.2500), SC_(19.625046917622010980803778160828770) }, { 1, SC_(-7.7500), SC_(19.618144334352289464741323510141514) }, { 1, SC_(-7.2500), SC_(19.610354539293269015698176691590937) }, { 1, SC_(-6.7500), SC_(19.601495010731061577124257953429755) }, { 1, SC_(-6.2500), SC_(19.591329569019785068016844943671793) }, { 1, SC_(-5.7500), SC_(19.579547136931335925546754524074474) }, { 1, SC_(-5.2500), SC_(19.565729569019785068016844943671793) }, { 1, SC_(-4.7500), SC_(19.549301390239464469970194977760674) }, { 1, SC_(-4.2500), SC_(19.529448389881463072551992335962043) }, { 1, SC_(-3.7500), SC_(19.504980060599575273294294700752364) }, { 1, SC_(-3.2500), SC_(19.474085068082155114074483685443011) }, { 1, SC_(-2.7500), SC_(19.433868949488464162183183589641253) }, { 1, SC_(-2.2500), SC_(19.379410511869137362595193744614609) }, { 1, SC_(-1.7500), SC_(19.301637544529786476232770366500757) }, { 1, SC_(-1.2500), SC_(19.181879647671606498397662880417078) }, { 1, SC_(-0.75000), SC_(18.975106932284888517049096897113002) }, { 1, SC_(-0.25000), SC_(18.541879647671606498397662880417078) },
{ 2, SC_(-12.750), SC_(-124.03079461415823384604153251543681) }, { 2, SC_(-12.250), SC_(124.01896466745858356132308878716344) }, { 2, SC_(-11.750), SC_(-124.03175955222881001960976796032603) }, { 2, SC_(-11.250), SC_(124.01787668541028735821044014586602) }, { 2, SC_(-10.750), SC_(-124.03299241970518808612682102178640) }, { 2, SC_(-10.250), SC_(124.01647202148710491650947992638728) }, { 2, SC_(-9.7500), SC_(-124.03460234084420729198290916496876) }, { 2, SC_(-9.2500), SC_(124.01461482266526541911391108670126) }, { 2, SC_(-8.7500), SC_(-124.03676016548723903560636876475972) }, { 2, SC_(-8.2500), SC_(124.01208782525148933477537240192445) }, { 2, SC_(-7.7500), SC_(-124.03974558822776381694747663647984) }, { 2, SC_(-7.2500), SC_(124.00852603656573370687098416695770) }, { 2, SC_(-6.7500), SC_(-124.04404218787195165891317097369578) }, { 2, SC_(-6.2500), SC_(124.00327776890408296268303058132483) }, { 2, SC_(-5.7500), SC_(-124.05054526159038888901020902683808) }, { 2, SC_(-5.2500), SC_(123.99508576890408296268303058132483) }, { 2, SC_(-4.7500), SC_(-124.06106552130930069964553408642549) }, { 2, SC_(-4.2500), SC_(123.98126436732757934536308673076874) }, { 2, SC_(-3.7500), SC_(-124.07972713378925404561433420306057) }, { 2, SC_(-3.2500), SC_(123.95521103942202265902072971875978) }, { 2, SC_(-2.7500), SC_(-124.11765305971517997154026012898649) }, { 2, SC_(-2.2500), SC_(123.89694977406016558118732052440384) }, { 2, SC_(-1.7500), SC_(-124.21382135423058192495874247308867) }, { 2, SC_(-1.2500), SC_(123.72136678366236036856729308956159) }, { 2, SC_(-0.75000), SC_(-124.58699919679617959259722643810325) }, { 2, SC_(-0.25000), SC_(122.69736678366236036856729308956159) },
@@ -158,7 +158,10 @@ void test_polygamma(T, const char* name)
{ SC_(3.0), SC_(-4.75), SC_(1558.5318783056006283387768251220856) }, { SC_(4.0), SC_(-4.75), SC_(-24481.582451925002662084791450344735) }, { SC_(5.0), SC_(-4.75), SC_(492231.26135104955223808370232711841) }, { SC_(6.0), SC_(-4.75), SC_(-1.1791224766632825018904254258839318e7) }, { SC_(7.0), SC_(-4.75), SC_(3.3035269584947086190323815702711063e8) }, { SC_(8.0), SC_(-4.75), SC_(-1.0569114259674567879524502974422220e10) }, { SC_(9.0), SC_(-4.75), SC_(3.8051374324232805216534792821594812e11) }, { SC_(10.0), SC_(-4.75), SC_(-1.5220204740668360644738342521168221e13) }, { SC_(11.0), SC_(-4.75), SC_(6.6969403856797200883878751610651220e14) }, { SC_(12.0), SC_(-4.75), SC_(-3.2145233093874118410246620594666049e16) }, { SC_(13.0), SC_(-4.75), SC_(1.6715535177261375554006502097069245e18) }, { SC_(14.0), SC_(-4.75), SC_(-9.3606970885707978200117071697879591e19) }, { SC_(15.0), SC_(-4.75), SC_(5.6164187748870728746282195405281448e21) }, { SC_(16.0), SC_(-4.75), SC_(-3.5945079045721164734032997263075953e23) }, { SC_(17.0), SC_(-4.75), SC_(2.4442654003427929640663758934254941e25) }, { SC_(18.0), SC_(-4.75), SC_(-1.7598710821897274459494051446575390e27) }, { SC_(19.0), SC_(-4.75), SC_(1.3375020239985042298043467743154361e29) }, { SC_(20.0), SC_(-4.75), SC_(-1.0700016187896297695358366297508349e31) }, { SC_(21.0), SC_(-4.75), SC_(8.9880135989785359476536358804280906e32) }, { SC_(22.0), SC_(-4.75), SC_(-7.9094519667650484338896180683559907e34) }, { SC_(23.0), SC_(-4.75), SC_(7.2766958095269026379022334905108714e36) },
{ SC_(1.0), SC_(-9.5), SC_(9.7696874450302318856305468284306792)}, { SC_(3.0), SC_(-9.5), SC_(194.81619198176773011863271713047162) }, { SC_(5.0), SC_(-9.5), SC_(15382.226860156995915624995579131219) }, { SC_(7.0), SC_(-9.5), SC_(2.5808804363008334969805587475917565e6) }, { SC_(9.0), SC_(-9.5), SC_(7.4319090477015599086877150154313919e8) }, { SC_(11.0), SC_(-9.5), SC_(3.2699904226951756570017589463296126e11) }, { SC_(13.0), SC_(-9.5), SC_(2.0404706026930390078103975418992502e14) }, { SC_(15.0), SC_(-9.5), SC_(1.7139949874533303450398622617829339e17) }, { SC_(17.0), SC_(-9.5), SC_(1.8648265078298896515398973583300966e20) }, { SC_(19.0), SC_(-9.5), SC_(2.5510826568574986072623191304028880e23) }, { SC_(21.0), SC_(-9.5), SC_(4.2858188624279670465725116159418790e26) },
{ 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_(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) }
{ 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");