Date: 2007-06-30 09:09:56
Date: 2007-06-30 09:09:55 EDT (Sat, 30 Jun 2007)
New Revision: 7321
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(-)
--- 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));
--- 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
--- 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 @@
-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
--- 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,
+"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
--- 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:
--- 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
+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,
--- 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
-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