diff options
| -rw-r--r-- | Doc/library/math.rst | 18 | ||||
| -rw-r--r-- | Lib/test/math_testcases.txt | 82 | ||||
| -rw-r--r-- | Misc/NEWS | 2 | ||||
| -rw-r--r-- | Modules/mathmodule.c | 142 |
4 files changed, 242 insertions, 2 deletions
diff --git a/Doc/library/math.rst b/Doc/library/math.rst index 39aea29..fce4d49 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -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 diff --git a/Lib/test/math_testcases.txt b/Lib/test/math_testcases.txt index bb99290..bd8813e 100644 --- a/Lib/test/math_testcases.txt +++ b/Lib/test/math_testcases.txt @@ -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. diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index ef84298..f835005 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -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}, |