
BoostCommit : 
From: johnmaddock_at_[hidden]
Date: 20070630 09:09:56
Author: johnmaddock
Date: 20070630 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 20070630 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 60term
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 20070630 09:09:55 EDT (Sat, 30 Jun 2007)
@@ 153,7 +153,7 @@
When the significand (mantissa) size is recognised
(currently for 53, 64 and 113bit reals, plus singleprecision 24bit 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 20070630 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 20070630 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 128bit long doubles for example), or because
+literature methods were unavailable or under nonBSL
+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 illconditioned (the coefficients can give a
+theoretically good solution, but the resulting rational function
+may be uncomputable 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 illdefined 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 20070630 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 128bit long doubles), rootpreserving
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 20070630 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.openstd.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.openstd.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 "realnumber" 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 20070630 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
BoostCommit 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