diff --git a/doc/digamma.qbk b/doc/digamma.qbk index 1ab430f6e..6d92d6fa9 100644 --- a/doc/digamma.qbk +++ b/doc/digamma.qbk @@ -40,9 +40,14 @@ than the one shown will have __zero_error. [[53] [Win32 Visual C++ 8] [Peak=0.98 Mean=0.36] [Peak=0.99 Mean=0.5] [Peak=0.95 Mean=0.5] [Peak=214 Mean=16] ] [[64] [Linux IA32 / GCC] [Peak=1.4 Mean=0.4] [Peak=1.3 Mean=0.45] [Peak=0.98 Mean=0.35] [Peak=180 Mean=13] ] [[64] [Linux IA64 / GCC] [Peak=0.92 Mean=0.4] [Peak=1.3 Mean=0.45] [Peak=0.98 Mean=0.4] [Peak=180 Mean=13] ] -[[113] [HPUX IA64, aCC A.06.06] [Peak=0.9 Mean=0.4] [Peak=200 Mean=14] [Peak=0.99 Mean=0.4] [Peak=200 Mean=14] ] +[[113] [HPUX IA64, aCC A.06.06] [Peak=0.9 Mean=0.4] [Peak=1.1 Mean=0.5] [Peak=0.99 Mean=0.4] [Peak=64 Mean=6] ] ] +As shown above, error rates for positive arguments are generally very low. +For negative arguments there are an infinite number of irrational roots: +relative errors very close to these can be arbitrarily large, although +absolute error will remain very low. + [h4 Testing] There are two sets of tests: spot values are computed using @@ -119,7 +124,7 @@ rational approximation. Note that since X0 is irrational, we need twice as many digits in X0 as in x in order to avoid cancellation error during the subtraction (this assumes that /x/ is an exact value, if it's not then all bets are off). That means that even when x is the value of the root rounded to the nearest -representable value, the result of digamma(x) ['[*will not be zero*]]. +representable value, the result of digamma(x) ['[*will not be zero]]. [endsect][/section:tgamma The Gamma Function] diff --git a/doc/equations/chf.png b/doc/equations/chf.png index e13dcf00d..2923d3137 100644 Binary files a/doc/equations/chf.png and b/doc/equations/chf.png differ diff --git a/doc/equations/dist_tutorial1.png b/doc/equations/dist_tutorial1.png index 80d54afe8..172eaab0b 100644 Binary files a/doc/equations/dist_tutorial1.png and b/doc/equations/dist_tutorial1.png differ diff --git a/doc/equations/dist_tutorial2.png b/doc/equations/dist_tutorial2.png index 8b2bf9348..3bf4ff048 100644 Binary files a/doc/equations/dist_tutorial2.png and b/doc/equations/dist_tutorial2.png differ diff --git a/doc/equations/dist_tutorial3.png b/doc/equations/dist_tutorial3.png index 205ba53a5..d6d5d8014 100644 Binary files a/doc/equations/dist_tutorial3.png and b/doc/equations/dist_tutorial3.png differ diff --git a/doc/equations/dist_tutorial4.png b/doc/equations/dist_tutorial4.png index 4a94b80a1..3e89ffb00 100644 Binary files a/doc/equations/dist_tutorial4.png and b/doc/equations/dist_tutorial4.png differ diff --git a/doc/equations/hazard.png b/doc/equations/hazard.png index 7a5978611..b65aebd46 100644 Binary files a/doc/equations/hazard.png and b/doc/equations/hazard.png differ diff --git a/doc/equations/students_t_dist.png b/doc/equations/students_t_dist.png index 5d4657d8c..3a0a239ae 100644 Binary files a/doc/equations/students_t_dist.png and b/doc/equations/students_t_dist.png differ diff --git a/doc/erf.qbk b/doc/erf.qbk index ccf337f81..82b4dda35 100644 --- a/doc/erf.qbk +++ b/doc/erf.qbk @@ -70,15 +70,6 @@ versions of these functions use differing implementations internally, so this gives us reasonably independent test data. Using our test data to test other "known good" implementations also provides an additional sanity check. -One should note that our tests rely on decimal to binary conversion of -floating-point numbers, and assume that the result will be correctly rounded. -In practice, it appears that in a few very rare cases the test data may be incorrect in the last bit: -this depends upon the compiler and standard library used, and means that the -relative errors quoted above have to treated somewhat circumspectly. Using -binary or hexadecimal coded test data would remove this issue (or at least -confirm whether it is actually an issue or not), but would make the test data -unportable, so is not used at present. - [h4 Implementation] All versions of these functions first use the usual reflection formulas diff --git a/include/boost/math/special_functions/digamma.hpp b/include/boost/math/special_functions/digamma.hpp index 0a925890e..a22f80c72 100644 --- a/include/boost/math/special_functions/digamma.hpp +++ b/include/boost/math/special_functions/digamma.hpp @@ -158,8 +158,8 @@ T digamma_imp_1_2(T x, const mpl::int_<0>*) static const T root1 = 1569415565.0 / 1073741824uL; static const T root2 = (381566830.0 / 1073741824uL) / 1073741824uL; static const T root3 = ((111616537.0 / 1073741824uL) / 1073741824uL) / 1073741824uL; - static const T root4 = (((503992063.0 / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL; - static const T root5 = 0.75231638452626400509999138382223723380394592598054e-36L; + static const T root4 = (((503992070.0 / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL; + static const T root5 = 0.52112228569249997894452490385577338504019838794544e-36L; static const T P[] = { 0.25479851061131551526977464225335883769L, diff --git a/minimax/f.cpp b/minimax/f.cpp index ed59131ce..862fd97c6 100644 --- a/minimax/f.cpp +++ b/minimax/f.cpp @@ -122,6 +122,10 @@ NTL::RR f(const NTL::RR& x, int variant) } case 10: { + // + // Digamma over the interval [1,2], set x-offset to 1 + // and optimise for absolute error over [0,1]. + // int current_precision = NTL::RR::precision(); if(current_precision < 1000) NTL::RR::SetPrecision(1000); @@ -129,11 +133,16 @@ NTL::RR f(const NTL::RR& x, int variant) // This value for the root of digamma is calculated using our // differentiated lanczos approximation. It agrees with Cody // to ~ 25 digits and to Morris to 35 digits. See: - // TOLMS ALGORITHM 708 (Didonato and Morris). + // TOMS ALGORITHM 708 (Didonato and Morris). // and Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher. // //NTL::RR root = boost::lexical_cast("1.4616321449683623412626595423257213234331845807102825466429633351908372838889871"); - + // + // Actually better to calculate the root on the fly, it appears to be more + // accurate: convergence is easier with the 1000-bit value, the approximation + // produced agrees with functions.mathworld.com values to 35 digits even quite + // near the root. + // static boost::math::tools::eps_tolerance tol(1000); static boost::uintmax_t max_iter = 1000; static const NTL::RR root = boost::math::tools::bracket_and_solve_root(&boost::math::digamma, NTL::RR(1.4), NTL::RR(1.5), true, tol, max_iter).first; @@ -142,6 +151,12 @@ NTL::RR f(const NTL::RR& x, int variant) double lim = 1e-65; if(fabs(x2 - root) < lim) { + // + // This is a problem area: + // x2-root suffers cancellation error, so does digamma. + // That gets compounded again when Remez calculates the error + // function. This cludge seems to stop the worst of the problems: + // static const NTL::RR a = boost::math::digamma(root - lim) / -lim; static const NTL::RR b = boost::math::digamma(root + lim) / lim; NTL::RR fract = (x2 - root + lim) / (2*lim); diff --git a/test/test_digamma.cpp b/test/test_digamma.cpp index 3f7dd4578..ecf244de0 100644 --- a/test/test_digamma.cpp +++ b/test/test_digamma.cpp @@ -42,14 +42,6 @@ void expected_results() // Define the max and mean errors expected for // various compilers and platforms. // - add_expected_result( - ".*aCC.*", // compiler - ".*", // stdlib - ".*", // platform - ".*", // test type(s) - ".*Root.*", // test data group - ".*", 220, 40); // test function - add_expected_result( ".*", // compiler ".*", // stdlib @@ -137,12 +129,11 @@ void test_spots(T, const char* t) // Special tolerance (200eps) for when we're very near the root, // and T has more than 64-bits in it's mantissa: // - T tolerance2 = boost::math::tools::digits() > 64 ? boost::math::tools::epsilon() * 20000 : tolerance; BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(0.125)), static_cast(-8.3884926632958548678027429230863430000514460424495L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(0.5)), static_cast(-1.9635100260214234794409763329987555671931596046604L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1)), static_cast(-0.57721566490153286060651209008240243104215933593992L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1.5)), static_cast(0.036489973978576520559023667001244432806840395339566L), tolerance); - BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1.5) - static_cast(1)/32), static_cast(0.00686541147073577672813890866512415766586241385896200579891429L), tolerance2); + BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(1.5) - static_cast(1)/32), static_cast(0.00686541147073577672813890866512415766586241385896200579891429L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(2)), static_cast(0.42278433509846713939348790991759756895784066406008L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(8)), static_cast(2.0156414779556099965363450527747404261006978069172L), tolerance); BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast(12)), static_cast(2.4426616799758120167383652547949424463027180089374L), tolerance);