libs/math/doc/sf/erf_inv.qbk

   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
  25 The return type of these functions is computed using the __arg_promotion_rules:
  26 the 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
  38 Returns the [@http://functions.wolfram.com/GammaBetaErf/InverseErf/ inverse error function]
  39 of 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
  51 Returns the inverse of the complement of the error function of z, that is a
  52 value x such that:
  53
  54 [expression ['p = erfc(x);]]
  55
  56 [graph erfc_inv]
  57
  58 [h4 Accuracy]
  59
  60 For types up to and including 80-bit long doubles the approximations used
  61 are accurate to less than ~ 2 epsilon.  For higher precision types these
  62 functions 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
  69 The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision,
  70 and 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
  80 There 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
  84 generous tolerances: in general both the error functions and their
  85 inverses change so rapidly in some places that round tripping to more than a couple
  86 of 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
  89 back where you started.
  90 * Accuracy checks using high-precision test values.  These measure
  91 the accuracy of the result, given /exact/ input values.
  92
  93 [h4 Implementation]
  94
  95 These functions use a rational approximation [jm_rationals]
  96 to calculate an initial
  97 approximation to the result that is accurate to ~10[super -19],
  98 then only if that has insufficient accuracy compared to the epsilon for T,
  99 do we clean up the result using
 100 [@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration].
 101
 102 Constructing rational approximations to the erf/erfc functions is actually
 103 surprisingly hard, especially at high precision.  For this reason no attempt
 104 has been made to achieve 10[super -34 ] accuracy suitable for use with 128-bit
 105 reals.
 106
 107 In the following discussion, /p/ is the value passed to erf_inv, and /q/ is
 108 the value passed to erfc_inv, so that /p = 1 - q/ and /q = 1 - p/ and in both
 109 cases we want to solve for the same result /x/.
 110
 111 For /p < 0.5/ the inverse erf function is reasonably smooth and the approximation:
 112
 113 [expression ['x = p(p + 10)(Y + R(p))]]
 114
 115 Gives a good result for a constant Y, and R(p) optimised for low absolute error
 116 compared to |Y|.
 117
 118 For q < 0.5 things get trickier, over the interval /0.5 > q > 0.25/
 119 the following approximation works well:
 120
 121 [expression ['x = sqrt(-2log(q)) / (Y + R(q))]]
 122
 123 While for q < 0.25, let
 124
 125 [expression ['z = sqrt(-log(q))]]
 126
 127 Then the result is given by:
 128
 129 [expression ['x = z(Y + R(z - B))]]
 130
 131 As before Y is a constant and the rational function R is optimised for low
 132 absolute error compared to |Y|.  B is also a constant: it is the smallest value
 133 of /z/ for which each approximation is valid.  There are several approximations
 134 of this form each of which reaches a little further into the tail of the erfc
 135 function (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 ]