1[section:error_inv Error Function Inverses] 2 3[h4 Synopsis] 4 5`` 6#include <boost/math/special_functions/erf.hpp> 7`` 8 9 namespace boost{ namespace math{ 10 11 template <class T> 12 ``__sf_result`` erf_inv(T p); 13 14 template <class T, class ``__Policy``> 15 ``__sf_result`` erf_inv(T p, const ``__Policy``&); 16 17 template <class T> 18 ``__sf_result`` erfc_inv(T p); 19 20 template <class T, class ``__Policy``> 21 ``__sf_result`` erfc_inv(T p, const ``__Policy``&); 22 23 }} // namespaces 24 25The return type of these functions is computed using the __arg_promotion_rules: 26the return type is `double` if T is an integer type, and T otherwise. 27 28[optional_policy] 29 30[h4 Description] 31 32 template <class T> 33 ``__sf_result`` erf_inv(T z); 34 35 template <class T, class ``__Policy``> 36 ``__sf_result`` erf_inv(T z, const ``__Policy``&); 37 38Returns the [@http://functions.wolfram.com/GammaBetaErf/InverseErf/ inverse error function] 39of z, that is a value x such that: 40 41[expression ['p = erf(x);]] 42 43[graph erf_inv] 44 45 template <class T> 46 ``__sf_result`` erfc_inv(T z); 47 48 template <class T, class ``__Policy``> 49 ``__sf_result`` erfc_inv(T z, const ``__Policy``&); 50 51Returns the inverse of the complement of the error function of z, that is a 52value x such that: 53 54[expression ['p = erfc(x);]] 55 56[graph erfc_inv] 57 58[h4 Accuracy] 59 60For types up to and including 80-bit long doubles the approximations used 61are accurate to less than ~ 2 epsilon. For higher precision types these 62functions have the same accuracy as the 63[link math_toolkit.sf_erf.error_function forward error functions]. 64 65[table_erf_inv] 66 67[table_erfc_inv] 68 69The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 70and GCC-7.1/Ubuntu for `long double` and `__float128`. 71 72[graph erfc__double] 73 74[graph erfc__80_bit_long_double] 75 76[graph erfc____float128] 77 78[h4 Testing] 79 80There are two sets of tests: 81 82* Basic sanity checks attempt to "round-trip" from 83/x/ to /p/ and back again. These tests have quite 84generous tolerances: in general both the error functions and their 85inverses change so rapidly in some places that round tripping to more than a couple 86of significant digits isn't possible. This is especially true when 87/p/ is very near one: in this case there isn't enough 88"information content" in the input to the inverse function to get 89back where you started. 90* Accuracy checks using high-precision test values. These measure 91the accuracy of the result, given /exact/ input values. 92 93[h4 Implementation] 94 95These functions use a rational approximation [jm_rationals] 96to calculate an initial 97approximation to the result that is accurate to ~10[super -19], 98then only if that has insufficient accuracy compared to the epsilon for T, 99do we clean up the result using 100[@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration]. 101 102Constructing rational approximations to the erf/erfc functions is actually 103surprisingly hard, especially at high precision. For this reason no attempt 104has been made to achieve 10[super -34 ] accuracy suitable for use with 128-bit 105reals. 106 107In the following discussion, /p/ is the value passed to erf_inv, and /q/ is 108the value passed to erfc_inv, so that /p = 1 - q/ and /q = 1 - p/ and in both 109cases we want to solve for the same result /x/. 110 111For /p < 0.5/ the inverse erf function is reasonably smooth and the approximation: 112 113[expression ['x = p(p + 10)(Y + R(p))]] 114 115Gives a good result for a constant Y, and R(p) optimised for low absolute error 116compared to |Y|. 117 118For q < 0.5 things get trickier, over the interval /0.5 > q > 0.25/ 119the following approximation works well: 120 121[expression ['x = sqrt(-2log(q)) / (Y + R(q))]] 122 123While for q < 0.25, let 124 125[expression ['z = sqrt(-log(q))]] 126 127Then the result is given by: 128 129[expression ['x = z(Y + R(z - B))]] 130 131As before Y is a constant and the rational function R is optimised for low 132absolute error compared to |Y|. B is also a constant: it is the smallest value 133of /z/ for which each approximation is valid. There are several approximations 134of this form each of which reaches a little further into the tail of the erfc 135function (at `long double` precision the extended exponent range compared to 136`double` means that the tail goes on for a very long way indeed). 137 138[endsect] [/ :error_inv The Error Function Inverses] 139 140[/ 141 Copyright 2006 John Maddock and Paul A. Bristow. 142 Distributed under the Boost Software License, Version 1.0. 143 (See accompanying file LICENSE_1_0.txt or copy at 144 http://www.boost.org/LICENSE_1_0.txt). 145] 146