(original) (raw)

Index: Misc/NEWS =================================================================== --- Misc/NEWS (revision 76762) +++ Misc/NEWS (working copy) @@ -1654,7 +1654,7 @@ - Issue #7078: Set struct.__doc__ from _struct.__doc__. -- Issue #3366: Add gamma, lgamma functions to math module. +- Issue #3366: Add erf, erfc, gamma, lgamma functions to math module. - Issue #6823: Allow time.strftime() to accept a tuple with a isdst field outside of the range of [-1, 1] by normalizing the value to within that Index: Doc/library/math.rst =================================================================== --- Doc/library/math.rst (revision 76762) +++ Doc/library/math.rst (working copy) @@ -311,6 +311,20 @@ Special functions ----------------- +.. function:: erf(x) + + Return the error function at *x*. + + .. versionadded:: 2.7 + + +.. function:: erfc(x) + + Return the complementary error function at *x*. + + .. versionadded:: 2.7 + + .. function:: gamma(x) Return the Gamma function at *x*. Index: Lib/erf.py =================================================================== --- Lib/erf.py (revision 0) +++ Lib/erf.py (revision 0) @@ -0,0 +1,114 @@ +# Implementations of the error and complementary error functions. + +from math import copysign, exp, isnan, isinf + +# 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; 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) < SERIES_CUTOFF and +# continued fraction expansion for SERIES_CUTOFF <= abs(x) < +# CONTFRAC_CUTOFF. SERIES_TERMS and CONTFRAC_TERMS are the default +# numbers of terms to use for the relevant expansions. +SERIES_CUTOFF = 1.5 +SERIES_TERMS = 25 +CONTFRAC_CUTOFF = 30.0 +CONTFRAC_TERMS = 50 + +SQRTPI = 1.7724538509055161 + +def erf_series(x, nterms=SERIES_TERMS): + """Error function, via power series. + + Given a finite float x, return an approximation to erf(x). + Converges reasonably fast for x small (say abs(x) < 1). + """ + + x2 = x*x + acc = 0.0 + fk = float(nterms) + 0.5 + for i in xrange(nterms): + acc = 2.0 + x2 * acc / fk + fk -= 1.0 + return x * exp(-x2) / SQRTPI * acc + +def erfc_cont_frac(x, nterms=CONTFRAC_TERMS): + """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. + + """ + + if x >= CONTFRAC_CUTOFF: + return 0.0 + + x2 = x * x + pkm1, qkm1 = 0.0, 1.0 + pk, qk = 1, 0.5 + x2 + for i in xrange(1, nterms): + ak = -i*(i-0.5) + bk = (4*i+1.0)/2.0 + x2 + pk, pkm1 = bk*pk + ak*pkm1, pk + qk, qkm1 = bk*qk + ak*qkm1, qk + return x * exp(-x2) / SQRTPI * (pk / qk) + +def erf(x): + """Error function of x.""" + + if isnan(x): + return x + + absx = abs(x) + if absx < SERIES_CUTOFF: + return erf_series(x) + else: + result = 1.0 - erfc_cont_frac(absx) + return result if x > 0.0 else -result + +def erfc(x): + """Complementary error function of x.""" + + if isnan(x): + return x + + absx = abs(x) + if absx < SERIES_CUTOFF: + return 1.0 - erf_series(x) + else: + result = erfc_cont_frac(absx) + return result if x > 0.0 else 2.0 - result + + +# Get corresponding correctly rounded functions from bigfloat module, +# for comparison. +#from bigfloat import erf as erf_bigfloat, erfc as erfc_bigfloat Index: Lib/test/test_math.py =================================================================== --- Lib/test/test_math.py (revision 76762) +++ Lib/test/test_math.py (working copy) @@ -4,6 +4,7 @@ from test.test_support import run_unittest, verbose import unittest import math +import erf import os import sys import random @@ -968,7 +969,10 @@ failures = [] for id, fn, arg, expected, flags in parse_mtestfile(math_testcases): - func = getattr(math, fn) + try: + func = getattr(math, fn) + except AttributeError: + func = getattr(erf, fn) if 'invalid' in flags or 'divide-by-zero' in flags: expected = 'ValueError' @@ -989,7 +993,7 @@ if not math.isnan(expected) and not math.isnan(got): # we use different closeness criteria for # different functions. - if fn == 'gamma': + if fn in ('gamma', 'erf', 'erfc'): accuracy_failure = ulps_check(expected, got, 20) elif fn == 'lgamma': accuracy_failure = acc_check(expected, got, Index: Lib/test/math_testcases.txt =================================================================== --- Lib/test/math_testcases.txt (revision 76762) +++ Lib/test/math_testcases.txt (working copy) @@ -47,6 +47,85 @@ -- 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 -- ---------------------------------------------------------