[Boost-commit] svn:boost r7321

From: johnmaddock_at_[hidden]
Date: 2007-06-30 09:09:56

Author: johnmaddock
Date: 2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
New Revision: 7321
URL: http://svn.boost.org/trac/boost/changeset/7321

Log:
Added rational approximation information.

Text files modified: 
   sandbox/math_toolkit/libs/math/doc/digamma.qbk        |     6 ++--                                    
   sandbox/math_toolkit/libs/math/doc/erf.qbk            |     2                                         
   sandbox/math_toolkit/libs/math/doc/erf_inv.qbk        |     3 +                                       
   sandbox/math_toolkit/libs/math/doc/implementation.qbk |    45 ++++++++++++++++++++++++++++++++++++--- 
   sandbox/math_toolkit/libs/math/doc/lgamma.qbk         |     2                                         
   sandbox/math_toolkit/libs/math/doc/math.qbk           |    15 ++++++++-----                           
   sandbox/math_toolkit/libs/math/doc/tgamma.qbk         |     2                                         
   7 files changed, 58 insertions(+), 17 deletions(-)

Modified: sandbox/math_toolkit/libs/math/doc/digamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/digamma.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/digamma.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -74,7 +74,7 @@
    
 is used to shift the evaluation to [1,2].
 
-For arguments in the range [1,2] a rational approximation is used (see below).
+For arguments in the range [1,2] a rational approximation [jm_rationals] is used (see below).
 
 For arguments in the range [2,BIG] the recurrence relation:
 
@@ -94,7 +94,7 @@
 the series to truncated after a suitably small number of terms and evaluated
 as a polynomial in `1/(x*x)`.
 
-The rational approximation in the range [1,2] is derived as follows.
+The rational approximation [jm_rationals] in the range [1,2] is derived as follows.
 
 First a high precision approximation to digamma was constructed using a 60-term
 differentiated __lanczos, the form used is:
@@ -118,7 +118,7 @@
 the value calculated by this method meets that requirement: the difficulty
 lies in independently verifying the value obtained.
 
-The rational approximation was optimised for absolute error using the form:
+The rational approximation [jm_rationals] was optimised for absolute error using the form:
 
    digamma(x) = (x - X0)(Y + R(x - 1));
    

Modified: sandbox/math_toolkit/libs/math/doc/erf.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/erf.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/erf.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -153,7 +153,7 @@
 
 When the significand (mantissa) size is recognised
 (currently for 53, 64 and 113-bit reals, plus single-precision 24-bit handled via promotion to double)
-then a series of rational approximations are used.
+then a series of rational approximations [jm_rationals] are used.
 
 For `z <= 0.5` then a rational approximation to erf is used, based on the 
 observation that:

Modified: sandbox/math_toolkit/libs/math/doc/erf_inv.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/erf_inv.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/erf_inv.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -67,7 +67,8 @@
 
 [h4 Implementation]
 
-These functions use a rational approximation to calculate an initial
+These functions use a rational approximation [jm_rationals] 
+to calculate an initial
 approximation to the result that is accurate to ~10[super -19], 
 then only if that has insufficient accuracy compared to the epsilon for T,
 do we clean up the result using

Modified: sandbox/math_toolkit/libs/math/doc/implementation.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/implementation.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/implementation.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -121,11 +121,9 @@
 *and* using try&catch is recommended for all applications.
 However, for simplicity, this is not done for most examples.]
 
-[h4 Handling of Functions that are Not Implemented]
+[h4 Handling of Functions that are Not Mathematically defined]
 
-Functions that are not implemented for any reason,
-usually because they are not defined, (or we were unable to find a definition),
-are handled as domain errors.
+Functions that are not mathematically defined are handled as domain errors.
 
 If the instruction 
 
@@ -215,6 +213,45 @@
 However as a approximation for the normal distribution,
 the most common usage, lower = -1, mode = 0 and upper = 1 would be more suitable.
 
+[h4 Rational Approximations Used]
+
+Some of the special functions in this library are implemented via
+rational approximations.  These are either taken from the literature,
+or devised by John Maddock using 
+[link math_toolkit.toolkit.internals2.minimax our Remez code].
+
+Rational rather than Polynomial approximations are used to ensure
+accuracy: polynomial approximations are often wonderful up to
+a certain level of accuracy, but then fail to provide much greater
+accuracy no matter how many more terms are added.
+
+Our own approximations were devised either for added accuracy
+(to support 128-bit long doubles for example), or because
+literature methods were unavailable or under non-BSL
+compatible license.  Our Remez code is known to produce good
+agreement with literature results in fairly simple "toy" cases,
+for more complex cases.  All approximations were checked 
+for convergence and to ensure that
+they were not ill-conditioned (the coefficients can give a 
+theoretically good solution, but the resulting rational function 
+may be un-computable at fixed precision).  
+
+Recomputing using different
+Remez implementations may well produce differing coefficients: the
+problem is well known to be ill conditioned in general, and our Remez implementation
+often found a broad and ill-defined minima for many of these approximations
+(of course for simple "toy" examples like approximating `exp` the minima 
+is well defined, and the coeffiecents should agree no matter whose Remez
+implementation is used).  This should not in general effect the validity
+of the approximations: there's good literature supporting the idea that 
+coefficients can be "in error" without necessarily adversely effecting 
+the result.  Note that "in error" has a special meaning in this context,
+see [@http://front.math.ucdavis.edu/0101.5042 
+"Approximate construction of rational approximations and the effect 
+of error autocorrection.", Grigori Litvinov, eprint arXiv:math/0101042].
+Therefore the coefficients still need to be accurately calculated, even if they can
+be in error compared to the "true" minimax solution.
+
 [h4 Representation of Mathematical Constants]
 
 A macro BOOST_DEFINE_MATH_CONSTANT in constants.hpp is used

Modified: sandbox/math_toolkit/libs/math/doc/lgamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/lgamma.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/lgamma.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -148,7 +148,7 @@
 
 For types with up to 113 bits of precision
 (up to and including 128-bit long doubles), root-preserving 
-rational approximations are used
+rational approximations [jm_rationals] are used
 over the intervals [1,2] and [2,3].  Over the interval [2,3] the approximation
 form used is:
 

Modified: sandbox/math_toolkit/libs/math/doc/math.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/math.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/math.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -45,6 +45,8 @@
                interfaces, library structure, and function and distribution names
                may be changed without notice.]
                
+[template tr1[] [@http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1836.pdf Technical Report on C++ Library Extensions]]
+[template jm_rationals[] [link math_toolkit.backgrounders.implementation.rational_approximations_used devised by JM]]
 [def __domain_error [link domain_error domain_error]]
 [def __pole_error [link pole_error pole_error]]
 [def __overflow_error [link overflow_error overflow_error]]
@@ -217,15 +219,16 @@
 Provides a small number of high quality 
 [@http://en.wikipedia.org/wiki/Special_functions special functions],
 initially these were concentrated on functions used in statistical applications
-but now also include some recently standardised in C99 & C++.
+along with those in the [tr1].
 
 The function families currently implemented are the gamma, beta & erf functions
 along with the incomplete gamma and beta functions (four variants
-of each) and all the possible inverses of these, plus digamma and various 
-special power and logarithmic functions.  Additions now include
-most new math functions added in the
-[@http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1836.pdf C++ Technical Report on Standard Library
-Extensions].
+of each) and all the possible inverses of these, plus digamma, 
+various factorial functions,
+Bessel functions, elliptic integrals, sinus cardinals (along with their 
+hyperbolic variants), inverse hyperbolic functions, Legrendre/Laguerre/Hermite
+polynomials and various 
+special power and logarithmic functions.
 
 All the implementations
 are fully generic and support the use of arbitrary "real-number" types,

Modified: sandbox/math_toolkit/libs/math/doc/tgamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/tgamma.qbk	(original)
+++ sandbox/math_toolkit/libs/math/doc/tgamma.qbk	2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
@@ -170,7 +170,7 @@
 Finally if the argument is a small integer then table lookup of the factorial
 is used.
 
-The function `tgamma1pm1` is implemented using rational approximations in the
+The function `tgamma1pm1` is implemented using rational approximations [jm_rationals] in the
 region `-0.5 < dz < 2`.  These are the same approximations (and internal routines)
 that are used for __lgamma, and so aren't detailed further here.  The result of
 the approximation is `log(tgamma(dz+1))` which can fed into __expm1 to give