Boost logo

Boost-Commit :

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


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk