]>
Commit | Line | Data |
---|---|---|
1 | [section:error_function Error Functions] | |
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(T z); | |
13 | ||
14 | template <class T, class ``__Policy``> | |
15 | ``__sf_result`` erf(T z, const ``__Policy``&); | |
16 | ||
17 | template <class T> | |
18 | ``__sf_result`` erfc(T z); | |
19 | ||
20 | template <class T, class ``__Policy``> | |
21 | ``__sf_result`` erfc(T z, 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(T z); | |
34 | ||
35 | template <class T, class ``__Policy``> | |
36 | ``__sf_result`` erf(T z, const ``__Policy``&); | |
37 | ||
38 | Returns the [@http://en.wikipedia.org/wiki/Error_function error function] | |
39 | [@http://functions.wolfram.com/GammaBetaErf/Erf/ erf] of z: | |
40 | ||
41 | [equation erf1] | |
42 | ||
43 | [graph erf] | |
44 | ||
45 | template <class T> | |
46 | ``__sf_result`` erfc(T z); | |
47 | ||
48 | template <class T, class ``__Policy``> | |
49 | ``__sf_result`` erfc(T z, const ``__Policy``&); | |
50 | ||
51 | Returns the complement of the [@http://functions.wolfram.com/GammaBetaErf/Erfc/ error function] of z: | |
52 | ||
53 | [equation erf2] | |
54 | ||
55 | [graph erfc] | |
56 | ||
57 | [h4 Accuracy] | |
58 | ||
59 | The following table shows the peak errors (in units of epsilon) | |
60 | found on various platforms with various floating point types, | |
61 | along with comparisons to the __gsl, __glibc, __hpc and __cephes libraries. | |
62 | Unless otherwise specified any floating point type that is narrower | |
63 | than the one shown will have __zero_error. | |
64 | ||
65 | [table_erf] | |
66 | ||
67 | [table_erfc] | |
68 | ||
69 | [h4 Testing] | |
70 | ||
71 | The tests for these functions come in two parts: | |
72 | basic sanity checks use spot values calculated using | |
73 | [@http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=Erf Mathworld's online evaluator], | |
74 | while accuracy checks use high-precision test values calculated at 1000-bit precision with | |
75 | [@http://shoup.net/ntl/doc/RR.txt NTL::RR] and this implementation. | |
76 | Note that the generic and type-specific | |
77 | versions of these functions use differing implementations internally, so this | |
78 | gives us reasonably independent test data. Using our test data to test other | |
79 | "known good" implementations also provides an additional sanity check. | |
80 | ||
81 | [h4 Implementation] | |
82 | ||
83 | All versions of these functions first use the usual reflection formulas | |
84 | to make their arguments positive: | |
85 | ||
86 | erf(-z) = 1 - erf(z); | |
87 | ||
88 | erfc(-z) = 2 - erfc(z); // preferred when -z < -0.5 | |
89 | ||
90 | erfc(-z) = 1 + erf(z); // preferred when -0.5 <= -z < 0 | |
91 | ||
92 | The generic versions of these functions are implemented in terms of | |
93 | the incomplete gamma function. | |
94 | ||
95 | When the significand (mantissa) size is recognised | |
96 | (currently for 53, 64 and 113-bit reals, plus single-precision 24-bit handled via promotion to double) | |
97 | then a series of rational approximations [jm_rationals] are used. | |
98 | ||
99 | For `z <= 0.5` then a rational approximation to erf is used, based on the | |
100 | observation that erf is an odd function and therefore erf is calculated using: | |
101 | ||
102 | erf(z) = z * (C + R(z*z)); | |
103 | ||
104 | where the rational approximation R(z*z) is optimised for absolute error: | |
105 | as long as its absolute error is small enough compared to the constant C, then any | |
106 | round-off error incurred during the computation of R(z*z) will effectively | |
107 | disappear from the result. As a result the error for erf and erfc in this | |
108 | region is very low: the last bit is incorrect in only a very small number of | |
109 | cases. | |
110 | ||
111 | For `z > 0.5` we observe that over a small interval \[a, b) then: | |
112 | ||
113 | erfc(z) * exp(z*z) * z ~ c | |
114 | ||
115 | for some constant c. | |
116 | ||
117 | Therefore for `z > 0.5` we calculate erfc using: | |
118 | ||
119 | erfc(z) = exp(-z*z) * (C + R(z - B)) / z; | |
120 | ||
121 | Again R(z - B) is optimised for absolute error, and the constant `C` is | |
122 | the average of `erfc(z) * exp(z*z) * z` taken at the endpoints of the range. | |
123 | Once again, as long as the absolute error in R(z - B) is small | |
124 | compared to `c` then `c + R(z - B)` will be correctly rounded, and the error | |
125 | in the result will depend only on the accuracy of the exp function. In practice, | |
126 | in all but a very small number of cases, the error is confined to the last bit | |
127 | of the result. The constant `B` is chosen so that the left hand end of the range | |
128 | of the rational approximation is 0. | |
129 | ||
130 | For large `z` over a range \[a, +[infin]\] the above approximation is modified to: | |
131 | ||
132 | erfc(z) = exp(-z*z) * (C + R(1 / z)) / z; | |
133 | ||
134 | [endsect] | |
135 | [/ :error_function The Error Functions] | |
136 | ||
137 | [/ | |
138 | Copyright 2006 John Maddock and Paul A. Bristow. | |
139 | Distributed under the Boost Software License, Version 1.0. | |
140 | (See accompanying file LICENSE_1_0.txt or copy at | |
141 | http://www.boost.org/LICENSE_1_0.txt). | |
142 | ] |