...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/special_functions/erf.hpp>
namespace boost{ namespace math{ template <class T> calculatedresulttype erf_inv(T p); template <class T, class Policy> calculatedresulttype erf_inv(T p, const Policy&); template <class T> calculatedresulttype erfc_inv(T p); template <class T, class Policy> calculatedresulttype erfc_inv(T p, const Policy&); }} // namespaces
The return type of these functions is computed using the result
type calculation rules: the return type is double
if T is an integer type, and T otherwise.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
template <class T> calculatedresulttype erf_inv(T z); template <class T, class Policy> calculatedresulttype erf_inv(T z, const Policy&);
Returns the inverse error function of z, that is a value x such that:
p = erf(x);
template <class T> calculatedresulttype erfc_inv(T z); template <class T, class Policy> calculatedresulttype 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);
For types up to and including 80bit long doubles the approximations used are accurate to less than ~ 2 epsilon. For higher precision types these functions have the same accuracy as the forward error functions.
Table 6.30. Error rates for erf_inv
Microsoft Visual C++ version 12.0 
GNU C++ version 5.1.0 
GNU C++ version 5.1.0 
Sun compiler version 0x5130 


Inverse Erf Function 
Max = 1.09ε (Mean = 0.502ε) 
Max = 0ε (Mean = 0ε) 
Max = 0.996ε (Mean = 0.389ε) 
Max = 0.996ε (Mean = 0.385ε) 
Table 6.31. Error rates for erfc_inv
Microsoft Visual C++ version 12.0 
GNU C++ version 5.1.0 
GNU C++ version 5.1.0 
Sun compiler version 0x5130 


Inverse Erfc Function 
Max = 1ε (Mean = 0.491ε) 
Max = 0ε (Mean = 0ε) 
Max = 0.996ε (Mean = 0.397ε) 
Max = 0.996ε (Mean = 0.397ε) 
Inverse Erfc Function: extreme values 
Max = 1.62ε (Mean = 0.383ε) 
Max = 1.62ε (Mean = 0.385ε) 
There are two sets of tests:
These functions use a rational approximation devised by JM to calculate an initial approximation to the result that is accurate to ~10^{19}, then only if that has insufficient accuracy compared to the epsilon for T, do we clean up the result using 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^{34 } accuracy suitable for use with 128bit 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).