[pypy-svn] r75951 - pypy/branch/fast-forward/pypy/module/math

[benjamin at codespeak.net]

Author: benjamin
Date: Wed Jul  7 05:39:42 2010
New Revision: 75951

Modified:
   pypy/branch/fast-forward/pypy/module/math/__init__.py
   pypy/branch/fast-forward/pypy/module/math/interp_math.py
Log:
add erf and erfc

Modified: pypy/branch/fast-forward/pypy/module/math/__init__.py
==============================================================================
--- pypy/branch/fast-forward/pypy/module/math/__init__.py	(original)
+++ pypy/branch/fast-forward/pypy/module/math/__init__.py	Wed Jul  7 05:39:42 2010
@@ -45,5 +45,7 @@
        'atanh'          : 'interp_math.atanh',
        'log1p'          : 'interp_math.log1p',
        'expm1'          : 'interp_math.expm1',
+       'erf'            : 'interp_math.erf',
+       'erfc'           : 'interp_math.erfc',
 }
 

Modified: pypy/branch/fast-forward/pypy/module/math/interp_math.py
==============================================================================
--- pypy/branch/fast-forward/pypy/module/math/interp_math.py	(original)
+++ pypy/branch/fast-forward/pypy/module/math/interp_math.py	Wed Jul  7 05:39:42 2010
@@ -421,3 +421,74 @@
     """exp(x) - 1"""
     return math1(space, rarithmetic.expm1, x)
 expm1.unwrap_spec = [ObjSpace, float]
+
+def erf(space, x):
+    """The error function"""
+    return math1(space, _erf, x)
+erf.unwrap_spec = [ObjSpace, float]
+
+def erfc(space, x):
+    """The complementary error function"""
+    return math1(space, _erfc, x)
+erfc.unwrap_spec = [ObjSpace, float]
+
+# Implementation of the error function, the complimentary error function, the
+# gamma function, and the natural log of the gamma function.  This exist in
+# libm, but I hear those implementations are horrible.
+
+ERF_SERIES_CUTOFF = 1.5
+ERF_SERIES_TERMS = 25
+ERFC_CONTFRAC_CUTOFF = 30.
+ERFC_CONTFRAC_TERMS = 50
+_sqrtpi = 1.772453850905516027298167483341145182798;
+
+def _erf_series(x):
+    x2 = x * x
+    acc = 0.
+    fk = ERF_SERIES_TERMS * + .5
+    for i in range(ERF_SERIES_TERMS):
+        acc = 2.0 * x2 * acc / fk
+        fk -= 1.
+    return acc * x * math.exp(-x2) / _sqrtpi
+
+def _erfc_contfrac(x):
+    if x >= ERFC_CONTFRAC_CUTOFF:
+        return 0.
+    x2 = x * x
+    a = 0.
+    da = 0.
+    p = 1.
+    p_last = 0.
+    q = da + x2
+    q_last = 1.
+    for i in range(ERFC_CONTFRAC_TERMS):
+        a += da
+        da += 2.
+        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 * math.exp(-x2) / _sqrtpi
+
+def _erf(x):
+    if rarithmetic.isnan(x):
+        return x
+    absx = abs(x)
+    if absx < ERF_SERIES_CUTOFF:
+        return _erf_series(x)
+    else:
+        cf = _erfc_contfrac(absx)
+        return 1. - cf if x > 0. else cf - 1.
+
+def _erfc(x):
+    if rarithmetic.isnan(x):
+        return x
+    absx = abs(x)
+    if absx < ERF_SERIES_CUTOFF:
+        return 1. - _erf_series(x)
+    else:
+        cf = _erfc_contfrac(absx)
+        return cf if x > 0. else 2. - cf