libs/math/doc/sf/erf.qbk

   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_pomotion_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 Errors In the Function erf(z)
  66 [[Significand Size] [Platform and Compiler] [z < 0.5]  [0.5 < z < 8] [z > 8]]
  67 [[53] [Win32, Visual C++ 8] [Peak=0 Mean=0
  68
  69 GSL Peak=2.0 Mean=0.3
  70
  71 __cephes Peak=1.1 Mean=0.7]  [Peak=0.9 Mean=0.09
  72
  73 GSL Peak=2.3 Mean=0.3
  74
  75 __cephes Peak=1.3 Mean=0.2] [Peak=0 Mean=0
  76
  77 GSL Peak=0 Mean=0
  78
  79 __cephes Peak=0 Mean=0]]
  80 [[64] [RedHat Linux IA32, gcc-3.3] [Peak=0.7 Mean=0.07
  81
  82 __glibc Peak=0.9 Mean=0.2]  [Peak=0.9 Mean=0.2
  83
  84 __glibc Peak=0.9 Mean=0.07] [Peak=0 Mean=0
  85
  86 __glibc Peak=0 Mean=0]]
  87 [[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=0.7 Mean=0.07
  88
  89 __glibc Peak=0 Mean=0]  [Peak=0.9 Mean=0.1
  90
  91 __glibc Peak=0.5 Mean=0.03] [Peak=0 Mean=0
  92
  93 __glibc Peak=0 Mean=0]]
  94 [[113] [HPUX IA64, aCC A.06.06] [Peak=0.8 Mean=0.1
  95
  96 __hpc Lib Peak=0.9 Mean=0.2]  [Peak=0.9 Mean=0.1
  97
  98 __hpc Lib Peak=0.5 Mean=0.02] [Peak=0 Mean=0
  99
 100 __hpc Lib Peak=0 Mean=0]]
 101 ]
 102
 103 [table Errors In the Function erfc(z)
 104 [[Significand Size] [Platform and Compiler] [z < 0.5]  [0.5 < z < 8] [z > 8]]
 105 [[53] [Win32, Visual C++ 8] [Peak=0.7 Mean=0.06
 106
 107 GSL Peak=1.0 Mean=0.4
 108
 109 __cephes Peak=0.7 Mean=0.06]  [Peak=0.99 Mean=0.3
 110
 111 GSL Peak=2.6 Mean=0.6
 112
 113 __cephes Peak=3.6 Mean=0.7] [Peak=1.0 Mean=0.2
 114
 115 GSL Peak=3.9 Mean=0.4
 116
 117 __cephes Peak=2.7 Mean=0.4]]
 118 [[64] [RedHat Linux IA32, gcc-3.3] [Peak=0 Mean=0
 119
 120 __glibc Peak=0 Mean=0]  [Peak=1.4 Mean=0.3
 121
 122 __glibc Peak=1.3 Mean=0.3] [Peak=1.6 Mean=0.4
 123
 124 __glibc Peak=1.3 Mean=0.4]]
 125 [[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=0 Mean=0
 126
 127 __glibc Peak=0 Mean=0]  [Peak=1.4 Mean=0.3
 128
 129 __glibc Peak=0 Mean=0] [Peak=1.5 Mean=0.4
 130
 131 __glibc Peak=0 Mean=0] ]
 132 [[113] [HPUX IA64, aCC A.06.06] [Peak=0 Mean=0
 133
 134 __hpc Peak=0 Mean=0]  [Peak=1.5 Mean=0.3
 135
 136 __hpc Peak=0.9 Mean=0.08] [Peak=1.6 Mean=0.4
 137
 138 __hpc Peak=0.9 Mean=0.1]]
 139 ]
 140
 141 [h4 Testing]
 142
 143 The tests for these functions come in two parts:
 144 basic sanity checks use spot values calculated using
 145 [@http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=Erf Mathworld's online evaluator],
 146 while accuracy checks use high-precision test values calculated at 1000-bit precision with
 147 [@http://shoup.net/ntl/doc/RR.txt NTL::RR] and this implementation.
 148 Note that the generic and type-specific
 149 versions of these functions use differing implementations internally, so this
 150 gives us reasonably independent test data.  Using our test data to test other
 151 "known good" implementations also provides an additional sanity check.
 152
 153 [h4 Implementation]
 154
 155 All versions of these functions first use the usual reflection formulas
 156 to make their arguments positive:
 157
 158    erf(-z) = 1 - erf(z);
 159
 160    erfc(-z) = 2 - erfc(z);  // preferred when -z < -0.5
 161
 162    erfc(-z) = 1 + erf(z);   // preferred when -0.5 <= -z < 0
 163
 164 The generic versions of these functions are implemented in terms of
 165 the incomplete gamma function.
 166
 167 When the significand (mantissa) size is recognised
 168 (currently for 53, 64 and 113-bit reals, plus single-precision 24-bit handled via promotion to double)
 169 then a series of rational approximations [jm_rationals] are used.
 170
 171 For `z <= 0.5` then a rational approximation to erf is used, based on the
 172 observation that erf is an odd function and therefore erf is calculated using:
 173
 174    erf(z) = z * (C + R(z*z));
 175
 176 where the rational approximation R(z*z) is optimised for absolute error:
 177 as long as its absolute error is small enough compared to the constant C, then any
 178 round-off error incurred during the computation of R(z*z) will effectively
 179 disappear from the result.  As a result the error for erf and erfc in this
 180 region is very low: the last bit is incorrect in only a very small number of
 181 cases.
 182
 183 For `z > 0.5` we observe that over a small interval \[a, b) then:
 184
 185    erfc(z) * exp(z*z) * z ~ c
 186
 187 for some constant c.
 188
 189 Therefore for `z > 0.5` we calculate erfc using:
 190
 191    erfc(z) = exp(-z*z) * (C + R(z - B)) / z;
 192
 193 Again R(z - B) is optimised for absolute error, and the constant `C` is
 194 the average of `erfc(z) * exp(z*z) * z` taken at the endpoints of the range.
 195 Once again, as long as the absolute error in R(z - B) is small
 196 compared to `c` then `c + R(z - B)` will be correctly rounded, and the error
 197 in the result will depend only on the accuracy of the exp function.  In practice,
 198 in all but a very small number of cases, the error is confined to the last bit
 199 of the result.  The constant `B` is chosen so that the left hand end of the range
 200 of the rational approximation is 0.
 201
 202 For large `z` over a range \[a, +[infin]\] the above approximation is modified to:
 203
 204    erfc(z) = exp(-z*z) * (C + R(1 / z)) / z;
 205
 206 [endsect]
 207 [/ :error_function The Error Functions]
 208
 209 [/
 210   Copyright 2006 John Maddock and Paul A. Bristow.
 211   Distributed under the Boost Software License, Version 1.0.
 212   (See accompanying file LICENSE_1_0.txt or copy at
 213   http://www.boost.org/LICENSE_1_0.txt).
 214 ]