123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150 |
- [section:error_function Error Function erf and complement erfc]
- [h4 Synopsis]
- ``
- #include <boost/math/special_functions/erf.hpp>
- ``
- namespace boost{ namespace math{
-
- template <class T>
- ``__sf_result`` erf(T z);
-
- template <class T, class ``__Policy``>
- ``__sf_result`` erf(T z, const ``__Policy``&);
-
- template <class T>
- ``__sf_result`` erfc(T z);
-
- template <class T, class ``__Policy``>
- ``__sf_result`` erfc(T z, const ``__Policy``&);
-
- }} // namespaces
-
- The return type of these functions is computed using the __arg_promotion_rules:
- the return type is `double` if T is an integer type, and T otherwise.
- [optional_policy]
- [h4 Description]
- template <class T>
- ``__sf_result`` erf(T z);
-
- template <class T, class ``__Policy``>
- ``__sf_result`` erf(T z, const ``__Policy``&);
-
- Returns the [@http://en.wikipedia.org/wiki/Error_function error function]
- [@http://functions.wolfram.com/GammaBetaErf/Erf/ erf] of z:
- [equation erf1]
- [graph erf]
- template <class T>
- ``__sf_result`` erfc(T z);
-
- template <class T, class ``__Policy``>
- ``__sf_result`` erfc(T z, const ``__Policy``&);
-
- Returns the complement of the [@http://functions.wolfram.com/GammaBetaErf/Erfc/ error function] of z:
- [equation erf2]
- [graph erfc]
- [h4 Accuracy]
- The following table shows the peak errors (in units of epsilon)
- found on various platforms with various floating-point types,
- along with comparisons to the __gsl, __glibc, __hpc and __cephes libraries.
- Unless otherwise specified any floating-point type that is narrower
- than the one shown will have __zero_error.
- [table_erf]
- [table_erfc]
- The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision,
- and GCC-7.1/Ubuntu for `long double` and `__float128`.
- [graph erf__double]
- [graph erf__80_bit_long_double]
- [graph erf____float128]
- [h4 Testing]
- The tests for these functions come in two parts:
- basic sanity checks use spot values calculated using
- [@http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=Erf Mathworld's online evaluator],
- while accuracy checks use high-precision test values calculated at 1000-bit precision with
- [@http://shoup.net/ntl/doc/RR.txt NTL::RR] and this implementation.
- Note that the generic and type-specific
- 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.
- [h4 Implementation]
- All versions of these functions first use the usual reflection formulas
- to make their arguments positive:
- [expression ['erf(-z) = 1 - erf(z);] ]
-
- [expression ['erfc(-z) = 2 - erfc(z);] // preferred when -z < -0.5]
-
- [expression ['erfc(-z) = 1 + erf(z);] // preferred when -0.5 <= -z < 0]
- The generic versions of these functions are implemented in terms of
- the incomplete gamma function.
- 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 [jm_rationals] are used.
- For `z <= 0.5` then a rational approximation to erf is used, based on the
- observation that erf is an odd function and therefore erf is calculated using:
- [expression ['erf(z) = z * (C + R(z*z));]]
-
- where the rational approximation /R(z*z)/ is optimised for absolute error:
- as long as its absolute error is small enough compared to the constant C, then any
- round-off error incurred during the computation of R(z*z) will effectively
- disappear from the result. As a result the error for erf and erfc in this
- region is very low: the last bit is incorrect in only a very small number of
- cases.
- For `z > 0.5` we observe that over a small interval \[['a, b)] then:
- [expression ['erfc(z) * exp(z*z) * z ~ c]]
-
- for some constant c.
- Therefore for `z > 0.5` we calculate `erfc` using:
- [expression ['erfc(z) = exp(-z*z) * (C + R(z - B)) / z;]]
-
- Again R(z - B) is optimised for absolute error, and the constant `C` is
- the average of `erfc(z) * exp(z*z) * z` taken at the endpoints of the range.
- Once again, as long as the absolute error in R(z - B) is small
- compared to `c` then `c + R(z - B)` will be correctly rounded, and the error
- in the result will depend only on the accuracy of the exp function. In practice,
- in all but a very small number of cases, the error is confined to the last bit
- of the result. The constant `B` is chosen so that the left hand end of the range
- of the rational approximation is 0.
- For large `z` over a range \[a, +[infin]\] the above approximation is modified to:
- [expression ['erfc(z) = exp(-z*z) * (C + R(1 / z)) / z;]]
- [endsect] [/section:error_function Error Function erf and complement erfc]
- [/
- Copyright 2006 John Maddock and Paul A. Bristow.
- Distributed under the Boost Software License, Version 1.0.
- (See accompanying file LICENSE_1_0.txt or copy at
- http://www.boost.org/LICENSE_1_0.txt).
- ]
|