]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | 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 | 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 | [h4 Testing] | |
70 | ||
71 | There are two sets of tests: | |
72 | ||
73 | * Basic sanity checks attempt to "round-trip" from | |
74 | /x/ to /p/ and back again. These tests have quite | |
75 | generous tolerances: in general both the error functions and their | |
76 | inverses change so rapidly in some places that round tripping to more than a couple | |
77 | of significant digits isn't possible. This is especially true when | |
78 | /p/ is very near one: in this case there isn't enough | |
79 | "information content" in the input to the inverse function to get | |
80 | back where you started. | |
81 | * Accuracy checks using high-precision test values. These measure | |
82 | the accuracy of the result, given /exact/ input values. | |
83 | ||
84 | [h4 Implementation] | |
85 | ||
86 | These functions use a rational approximation [jm_rationals] | |
87 | to calculate an initial | |
88 | approximation to the result that is accurate to ~10[super -19], | |
89 | then only if that has insufficient accuracy compared to the epsilon for T, | |
90 | do we clean up the result using | |
91 | [@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration]. | |
92 | ||
93 | Constructing rational approximations to the erf/erfc functions is actually | |
94 | surprisingly hard, especially at high precision. For this reason no attempt | |
95 | has been made to achieve 10[super -34 ] accuracy suitable for use with 128-bit | |
96 | reals. | |
97 | ||
98 | In the following discussion, /p/ is the value passed to erf_inv, and /q/ is | |
99 | the value passed to erfc_inv, so that /p = 1 - q/ and /q = 1 - p/ and in both | |
100 | cases we want to solve for the same result /x/. | |
101 | ||
102 | For /p < 0.5/ the inverse erf function is reasonably smooth and the approximation: | |
103 | ||
104 | x = p(p + 10)(Y + R(p)) | |
105 | ||
106 | Gives a good result for a constant Y, and R(p) optimised for low absolute error | |
107 | compared to |Y|. | |
108 | ||
109 | For q < 0.5 things get trickier, over the interval /0.5 > q > 0.25/ | |
110 | the following approximation works well: | |
111 | ||
112 | x = sqrt(-2log(q)) / (Y + R(q)) | |
113 | ||
114 | While for q < 0.25, let | |
115 | ||
116 | z = sqrt(-log(q)) | |
117 | ||
118 | Then the result is given by: | |
119 | ||
120 | x = z(Y + R(z - B)) | |
121 | ||
122 | As before Y is a constant and the rational function R is optimised for low | |
123 | absolute error compared to |Y|. B is also a constant: it is the smallest value | |
124 | of /z/ for which each approximation is valid. There are several approximations | |
125 | of this form each of which reaches a little further into the tail of the erfc | |
126 | function (at `long double` precision the extended exponent range compared to | |
127 | `double` means that the tail goes on for a very long way indeed). | |
128 | ||
129 | [endsect][/ :error_inv The Error Function Inverses] | |
130 | ||
131 | [/ | |
132 | Copyright 2006 John Maddock and Paul A. Bristow. | |
133 | Distributed under the Boost Software License, Version 1.0. | |
134 | (See accompanying file LICENSE_1_0.txt or copy at | |
135 | http://www.boost.org/LICENSE_1_0.txt). | |
136 | ] |