From e70f53e670a68ac704ea2382261e855bbe7ca29d Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 31 Oct 2014 18:51:14 +0000 Subject: [PATCH] [Polygamma] Fix issues with small negative arguments for x. Tidy up a couple of cosmetic coding issues. --- .../math/special_functions/detail/polygamma.hpp | 13 ++++++------- test/test_polygamma.hpp | 7 +++++-- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/include/boost/math/special_functions/detail/polygamma.hpp b/include/boost/math/special_functions/detail/polygamma.hpp index 5b94237fd..50c581686 100644 --- a/include/boost/math/special_functions/detail/polygamma.hpp +++ b/include/boost/math/special_functions/detail/polygamma.hpp @@ -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()); - part_term += log(n * (n + 1)) - boost::math::constants::ln_two() - log(x); + part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two() - 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(0.4F * policies::digits_base10()); - const int N4dn = static_cast(d4d + (4 * n)); - const int N = static_cast((std::min)(N4dn, (std::numeric_limits::max)())); + const int d4d = static_cast(0.4F * policies::digits_base10()); + const int N = d4d + (4 * n); const int m = n; const int iter = N - itrunc(x); @@ -234,7 +233,7 @@ } template - 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() * poly_cot_pi(n, z, pol, function); + T result = polygamma_imp(n, z, pol) + constants::pi() * poly_cot_pi(n, z, x, pol, function); return n & 1 ? T(-result) : result; } // diff --git a/test/test_polygamma.hpp b/test/test_polygamma.hpp index 335f715d0..da0cfa153 100644 --- a/test/test_polygamma.hpp +++ b/test/test_polygamma.hpp @@ -137,7 +137,7 @@ void test_polygamma(T, const char* name) } }; do_test_polygamma(big_data, name, "Mathematica Data - large arguments"); - boost::array, 373> neg_data = + boost::array, 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(neg_data, name, "Mathematica Data - negative arguments");