Fixed digamma 128-bit error rates with a more accurate value for the root.

Updated png's: hopefully they'll all render OK now.
Tweeked docs to reflect recent advances.


[SVN r3237]

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]

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

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,

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<NTL::RR>("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<NTL::RR> 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);

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<T>() > 64 ? boost::math::tools::epsilon<T>() * 20000 : tolerance;
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(0.125)), static_cast<T>(-8.3884926632958548678027429230863430000514460424495L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(0.5)), static_cast<T>(-1.9635100260214234794409763329987555671931596046604L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(1)), static_cast<T>(-0.57721566490153286060651209008240243104215933593992L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(1.5)), static_cast<T>(0.036489973978576520559023667001244432806840395339566L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(1.5) - static_cast<T>(1)/32), static_cast<T>(0.00686541147073577672813890866512415766586241385896200579891429L), tolerance2);
+   BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(1.5) - static_cast<T>(1)/32), static_cast<T>(0.00686541147073577672813890866512415766586241385896200579891429L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(2)), static_cast<T>(0.42278433509846713939348790991759756895784066406008L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(8)), static_cast<T>(2.0156414779556099965363450527747404261006978069172L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::digamma(static_cast<T>(12)), static_cast<T>(2.4426616799758120167383652547949424463027180089374L), tolerance);