| [section:error_inv Error Function Inverses] |
| |
| [h4 Synopsis] |
| |
| `` |
| #include <boost/math/special_functions/erf.hpp> |
| `` |
| |
| namespace boost{ namespace math{ |
| |
| template <class T> |
| ``__sf_result`` erf_inv(T p); |
| |
| template <class T, class ``__Policy``> |
| ``__sf_result`` erf_inv(T p, const ``__Policy``&); |
| |
| template <class T> |
| ``__sf_result`` erfc_inv(T p); |
| |
| template <class T, class ``__Policy``> |
| ``__sf_result`` erfc_inv(T p, const ``__Policy``&); |
| |
| }} // namespaces |
| |
| The return type of these functions is computed using the __arg_pomotion_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_inv(T z); |
| |
| template <class T, class ``__Policy``> |
| ``__sf_result`` erf_inv(T z, const ``__Policy``&); |
| |
| Returns the [@http://functions.wolfram.com/GammaBetaErf/InverseErf/ inverse error function] |
| of z, that is a value x such that: |
| |
| p = erf(x); |
| |
| [graph erf_inv] |
| |
| template <class T> |
| ``__sf_result`` erfc_inv(T z); |
| |
| template <class T, class ``__Policy``> |
| ``__sf_result`` erfc_inv(T z, const ``__Policy``&); |
| |
| Returns the inverse of the complement of the error function of z, that is a |
| value x such that: |
| |
| p = erfc(x); |
| |
| [graph erfc_inv] |
| |
| [h4 Accuracy] |
| |
| For types up to and including 80-bit long doubles the approximations used |
| are accurate to less than ~ 2 epsilon. For higher precision types these |
| functions have the same accuracy as the |
| [link math_toolkit.special.sf_erf.error_function forward error functions]. |
| |
| [h4 Testing] |
| |
| There are two sets of tests: |
| |
| * Basic sanity checks attempt to "round-trip" from |
| /x/ to /p/ and back again. These tests have quite |
| generous tolerances: in general both the error functions and their |
| inverses change so rapidly in some places that round tripping to more than a couple |
| of significant digits isn't possible. This is especially true when |
| /p/ is very near one: in this case there isn't enough |
| "information content" in the input to the inverse function to get |
| back where you started. |
| * Accuracy checks using high-precision test values. These measure |
| the accuracy of the result, given /exact/ input values. |
| |
| [h4 Implementation] |
| |
| 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 |
| [@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration]. |
| |
| Constructing rational approximations to the erf/erfc functions is actually |
| surprisingly hard, especially at high precision. For this reason no attempt |
| has been made to achieve 10[super -34 ] accuracy suitable for use with 128-bit |
| reals. |
| |
| In the following discussion, /p/ is the value passed to erf_inv, and /q/ is |
| the value passed to erfc_inv, so that /p = 1 - q/ and /q = 1 - p/ and in both |
| cases we want to solve for the same result /x/. |
| |
| For /p < 0.5/ the inverse erf function is reasonably smooth and the approximation: |
| |
| x = p(p + 10)(Y + R(p)) |
| |
| Gives a good result for a constant Y, and R(p) optimised for low absolute error |
| compared to |Y|. |
| |
| For q < 0.5 things get trickier, over the interval /0.5 > q > 0.25/ |
| the following approximation works well: |
| |
| x = sqrt(-2log(q)) / (Y + R(q)) |
| |
| While for q < 0.25, let |
| |
| z = sqrt(-log(q)) |
| |
| Then the result is given by: |
| |
| x = z(Y + R(z - B)) |
| |
| As before Y is a constant and the rational function R is optimised for low |
| absolute error compared to |Y|. B is also a constant: it is the smallest value |
| of /z/ for which each approximation is valid. There are several approximations |
| of this form each of which reaches a little further into the tail of the erfc |
| function (at `long double` precision the extended exponent range compared to |
| `double` means that the tail goes on for a very long way indeed). |
| |
| [endsect][/ :error_inv The Error Function Inverses] |
| |
| [/ |
| 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). |
| ] |