Merged revisions 76878 via svnmerge from
svn+ssh://pythondev@svn.python.org/python/trunk ........ r76878 | mark.dickinson | 2009-12-19 11:07:23 +0000 (Sat, 19 Dec 2009) | 3 lines Issue #3366: Add error function and complementary error function to math module. ........
This commit is contained in:
Mark Dickinson
parent
056aafe7f2
commit
45f992a45e
|
|
@ -161,6 +161,8 @@ Power and logarithmic functions
|
|||
>>> expm1(1e-5) # result accurate to full precision
|
||||
1.0000050000166668e-05
|
||||
|
||||
.. versionadded:: 3.2
|
||||
|
||||
|
||||
.. function:: log(x[, base])
|
||||
|
||||
|
|
@ -295,6 +297,20 @@ Hyperbolic functions
|
|||
Special functions
|
||||
-----------------
|
||||
|
||||
.. function:: erf(x)
|
||||
|
||||
Return the error function at *x*.
|
||||
|
||||
.. versionadded:: 3.2
|
||||
|
||||
|
||||
.. function:: erfc(x)
|
||||
|
||||
Return the complementary error function at *x*.
|
||||
|
||||
.. versionadded:: 3.2
|
||||
|
||||
|
||||
.. function:: gamma(x)
|
||||
|
||||
Return the Gamma function at *x*.
|
||||
|
|
@ -307,7 +323,7 @@ Special functions
|
|||
Return the natural logarithm of the absolute value of the Gamma
|
||||
function at *x*.
|
||||
|
||||
.. versionadded:: 2.7
|
||||
.. versionadded:: 3.2
|
||||
|
||||
|
||||
Constants
|
||||
|
|
|
|||
|
|
@ -47,6 +47,87 @@
|
|||
-- MPFR homepage at http://www.mpfr.org for more information about the
|
||||
-- MPFR project.
|
||||
|
||||
|
||||
-------------------------
|
||||
-- erf: error function --
|
||||
-------------------------
|
||||
|
||||
erf0000 erf 0.0 -> 0.0
|
||||
erf0001 erf -0.0 -> -0.0
|
||||
erf0002 erf inf -> 1.0
|
||||
erf0003 erf -inf -> -1.0
|
||||
erf0004 erf nan -> nan
|
||||
|
||||
-- tiny values
|
||||
erf0010 erf 1e-308 -> 1.1283791670955125e-308
|
||||
erf0011 erf 5e-324 -> 4.9406564584124654e-324
|
||||
erf0012 erf 1e-10 -> 1.1283791670955126e-10
|
||||
|
||||
-- small integers
|
||||
erf0020 erf 1 -> 0.84270079294971489
|
||||
erf0021 erf 2 -> 0.99532226501895271
|
||||
erf0022 erf 3 -> 0.99997790950300136
|
||||
erf0023 erf 4 -> 0.99999998458274209
|
||||
erf0024 erf 5 -> 0.99999999999846256
|
||||
erf0025 erf 6 -> 1.0
|
||||
|
||||
erf0030 erf -1 -> -0.84270079294971489
|
||||
erf0031 erf -2 -> -0.99532226501895271
|
||||
erf0032 erf -3 -> -0.99997790950300136
|
||||
erf0033 erf -4 -> -0.99999998458274209
|
||||
erf0034 erf -5 -> -0.99999999999846256
|
||||
erf0035 erf -6 -> -1.0
|
||||
|
||||
-- huge values should all go to +/-1, depending on sign
|
||||
erf0040 erf -40 -> -1.0
|
||||
erf0041 erf 1e16 -> 1.0
|
||||
erf0042 erf -1e150 -> -1.0
|
||||
erf0043 erf 1.7e308 -> 1.0
|
||||
|
||||
|
||||
----------------------------------------
|
||||
-- erfc: complementary error function --
|
||||
----------------------------------------
|
||||
|
||||
erfc0000 erfc 0.0 -> 1.0
|
||||
erfc0001 erfc -0.0 -> 1.0
|
||||
erfc0002 erfc inf -> 0.0
|
||||
erfc0003 erfc -inf -> 2.0
|
||||
erfc0004 erfc nan -> nan
|
||||
|
||||
-- tiny values
|
||||
erfc0010 erfc 1e-308 -> 1.0
|
||||
erfc0011 erfc 5e-324 -> 1.0
|
||||
erfc0012 erfc 1e-10 -> 0.99999999988716204
|
||||
|
||||
-- small integers
|
||||
erfc0020 erfc 1 -> 0.15729920705028513
|
||||
erfc0021 erfc 2 -> 0.0046777349810472662
|
||||
erfc0022 erfc 3 -> 2.2090496998585441e-05
|
||||
erfc0023 erfc 4 -> 1.541725790028002e-08
|
||||
erfc0024 erfc 5 -> 1.5374597944280349e-12
|
||||
erfc0025 erfc 6 -> 2.1519736712498913e-17
|
||||
|
||||
erfc0030 erfc -1 -> 1.8427007929497148
|
||||
erfc0031 erfc -2 -> 1.9953222650189528
|
||||
erfc0032 erfc -3 -> 1.9999779095030015
|
||||
erfc0033 erfc -4 -> 1.9999999845827421
|
||||
erfc0034 erfc -5 -> 1.9999999999984626
|
||||
erfc0035 erfc -6 -> 2.0
|
||||
|
||||
-- as x -> infinity, erfc(x) behaves like exp(-x*x)/x/sqrt(pi)
|
||||
erfc0040 erfc 20 -> 5.3958656116079012e-176
|
||||
erfc0041 erfc 25 -> 8.3001725711965228e-274
|
||||
erfc0042 erfc 27 -> 5.2370464393526292e-319
|
||||
erfc0043 erfc 28 -> 0.0
|
||||
|
||||
-- huge values
|
||||
erfc0050 erfc -40 -> 2.0
|
||||
erfc0051 erfc 1e16 -> 0.0
|
||||
erfc0052 erfc -1e150 -> 2.0
|
||||
erfc0053 erfc 1.7e308 -> 0.0
|
||||
|
||||
|
||||
---------------------------------------------------------
|
||||
-- lgamma: log of absolute value of the gamma function --
|
||||
---------------------------------------------------------
|
||||
|
|
@ -250,6 +331,7 @@ gam0132 gamma -4503599627370495.5 -> 0.0
|
|||
gam0140 gamma -63.349078729022985 -> 4.1777971677761880e-88
|
||||
gam0141 gamma -127.45117632943295 -> 1.1831110896236810e-214
|
||||
|
||||
|
||||
-----------------------------------------------------------
|
||||
-- expm1: exp(x) - 1, without precision loss for small x --
|
||||
-----------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -457,7 +457,7 @@ Extension Modules
|
|||
|
||||
- Issue #7078: Set struct.__doc__ from _struct.__doc__.
|
||||
|
||||
- Issue #3366: Add expm1, gamma, lgamma functions to math module.
|
||||
- Issue #3366: Add erf, erfc, expm1, gamma, lgamma functions to math module.
|
||||
|
||||
- Issue #6877: It is now possible to link the readline extension to the
|
||||
libedit readline emulation on OSX 10.5 or later.
|
||||
|
|
|
|||
|
|
@ -69,6 +69,7 @@ extern double copysign(double, double);
|
|||
*/
|
||||
|
||||
static const double pi = 3.141592653589793238462643383279502884197;
|
||||
static const double sqrtpi = 1.772453850905516027298167483341145182798;
|
||||
|
||||
static double
|
||||
sinpi(double x)
|
||||
|
|
@ -375,6 +376,141 @@ m_lgamma(double x)
|
|||
return r;
|
||||
}
|
||||
|
||||
/*
|
||||
Implementations of the error function erf(x) and the complementary error
|
||||
function erfc(x).
|
||||
|
||||
Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed.,
|
||||
Cambridge University Press), we use a series approximation for erf for
|
||||
small x, and a continued fraction approximation for erfc(x) for larger x;
|
||||
combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
|
||||
this gives us erf(x) and erfc(x) for all x.
|
||||
|
||||
The series expansion used is:
|
||||
|
||||
erf(x) = x*exp(-x*x)/sqrt(pi) * [
|
||||
2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
|
||||
|
||||
The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
|
||||
This series converges well for smallish x, but slowly for larger x.
|
||||
|
||||
The continued fraction expansion used is:
|
||||
|
||||
erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
|
||||
3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
|
||||
|
||||
after the first term, the general term has the form:
|
||||
|
||||
k*(k-0.5)/(2*k+0.5 + x**2 - ...).
|
||||
|
||||
This expansion converges fast for larger x, but convergence becomes
|
||||
infinitely slow as x approaches 0.0. The (somewhat naive) continued
|
||||
fraction evaluation algorithm used below also risks overflow for large x;
|
||||
but for large x, erfc(x) == 0.0 to within machine precision. (For
|
||||
example, erfc(30.0) is approximately 2.56e-393).
|
||||
|
||||
Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
|
||||
continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
|
||||
ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
|
||||
numbers of terms to use for the relevant expansions. */
|
||||
|
||||
#define ERF_SERIES_CUTOFF 1.5
|
||||
#define ERF_SERIES_TERMS 25
|
||||
#define ERFC_CONTFRAC_CUTOFF 30.0
|
||||
#define ERFC_CONTFRAC_TERMS 50
|
||||
|
||||
/*
|
||||
Error function, via power series.
|
||||
|
||||
Given a finite float x, return an approximation to erf(x).
|
||||
Converges reasonably fast for small x.
|
||||
*/
|
||||
|
||||
static double
|
||||
m_erf_series(double x)
|
||||
{
|
||||
double x2, acc, fk;
|
||||
int i;
|
||||
|
||||
x2 = x * x;
|
||||
acc = 0.0;
|
||||
fk = (double)ERF_SERIES_TERMS + 0.5;
|
||||
for (i = 0; i < ERF_SERIES_TERMS; i++) {
|
||||
acc = 2.0 + x2 * acc / fk;
|
||||
fk -= 1.0;
|
||||
}
|
||||
return acc * x * exp(-x2) / sqrtpi;
|
||||
}
|
||||
|
||||
/*
|
||||
Complementary error function, via continued fraction expansion.
|
||||
|
||||
Given a positive float x, return an approximation to erfc(x). Converges
|
||||
reasonably fast for x large (say, x > 2.0), and should be safe from
|
||||
overflow if x and nterms are not too large. On an IEEE 754 machine, with x
|
||||
<= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
|
||||
than the smallest representable nonzero float. */
|
||||
|
||||
static double
|
||||
m_erfc_contfrac(double x)
|
||||
{
|
||||
double x2, a, da, p, p_last, q, q_last, b;
|
||||
int i;
|
||||
|
||||
if (x >= ERFC_CONTFRAC_CUTOFF)
|
||||
return 0.0;
|
||||
|
||||
x2 = x*x;
|
||||
a = 0.0;
|
||||
da = 0.5;
|
||||
p = 1.0; p_last = 0.0;
|
||||
q = da + x2; q_last = 1.0;
|
||||
for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
|
||||
double temp;
|
||||
a += da;
|
||||
da += 2.0;
|
||||
b = da + x2;
|
||||
temp = p; p = b*p - a*p_last; p_last = temp;
|
||||
temp = q; q = b*q - a*q_last; q_last = temp;
|
||||
}
|
||||
return p / q * x * exp(-x2) / sqrtpi;
|
||||
}
|
||||
|
||||
/* Error function erf(x), for general x */
|
||||
|
||||
static double
|
||||
m_erf(double x)
|
||||
{
|
||||
double absx, cf;
|
||||
|
||||
if (Py_IS_NAN(x))
|
||||
return x;
|
||||
absx = fabs(x);
|
||||
if (absx < ERF_SERIES_CUTOFF)
|
||||
return m_erf_series(x);
|
||||
else {
|
||||
cf = m_erfc_contfrac(absx);
|
||||
return x > 0.0 ? 1.0 - cf : cf - 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
/* Complementary error function erfc(x), for general x. */
|
||||
|
||||
static double
|
||||
m_erfc(double x)
|
||||
{
|
||||
double absx, cf;
|
||||
|
||||
if (Py_IS_NAN(x))
|
||||
return x;
|
||||
absx = fabs(x);
|
||||
if (absx < ERF_SERIES_CUTOFF)
|
||||
return 1.0 - m_erf_series(x);
|
||||
else {
|
||||
cf = m_erfc_contfrac(absx);
|
||||
return x > 0.0 ? cf : 2.0 - cf;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
wrapper for atan2 that deals directly with special cases before
|
||||
|
|
@ -721,6 +857,10 @@ FUNC1(cos, cos, 0,
|
|||
"cos(x)\n\nReturn the cosine of x (measured in radians).")
|
||||
FUNC1(cosh, cosh, 1,
|
||||
"cosh(x)\n\nReturn the hyperbolic cosine of x.")
|
||||
FUNC1A(erf, m_erf,
|
||||
"erf(x)\n\nError function at x.")
|
||||
FUNC1A(erfc, m_erfc,
|
||||
"erfc(x)\n\nComplementary error function at x.")
|
||||
FUNC1(exp, exp, 1,
|
||||
"exp(x)\n\nReturn e raised to the power of x.")
|
||||
FUNC1(expm1, m_expm1, 1,
|
||||
|
|
@ -1497,6 +1637,8 @@ static PyMethodDef math_methods[] = {
|
|||
{"cos", math_cos, METH_O, math_cos_doc},
|
||||
{"cosh", math_cosh, METH_O, math_cosh_doc},
|
||||
{"degrees", math_degrees, METH_O, math_degrees_doc},
|
||||
{"erf", math_erf, METH_O, math_erf_doc},
|
||||
{"erfc", math_erfc, METH_O, math_erfc_doc},
|
||||
{"exp", math_exp, METH_O, math_exp_doc},
|
||||
{"expm1", math_expm1, METH_O, math_expm1_doc},
|
||||
{"fabs", math_fabs, METH_O, math_fabs_doc},
|
||||
|
|
|
|||
Loading…
Reference in New Issue