# Ticket #1173: trac_1173_complex_erf_v3.patch

File trac_1173_complex_erf_v3.patch, 7.4 KB (added by dsm, 11 years ago)

11143-style erf, use mpfr where possible

• ## sage/functions/other.py

# HG changeset patch
# User D. S. McNeil <dsm054@gmail.com>
# Date 1320009980 14400
# Node ID f123a88b4f42992eff416e8af743bbb2a3f73fea
# Parent  3b825be5f6ecb07d743c9b76c75fabd0df670740
Trac #1173: add support for complex erf

diff --git a/sage/functions/other.py b/sage/functions/other.py
 a from sage.symbolic.pynac import py_factorial_py from sage.libs.pari.gen import pari from sage.symbolic.all import SR from sage.rings.all import Integer, Rational, RealField, RR, \ from sage.rings.all import Integer, Rational, RealField, RR, RDF, \ is_ComplexNumber, ComplexField from sage.misc.latex import latex import math one_half = ~SR(2) class Function_erf(BuiltinFunction): _eval_ = BuiltinFunction._eval_default r""" The error function, defined for real values as \text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt. This function is also defined for complex values, via analytic continuation. EXAMPLES: We can evaluate numerically via mpmath:: sage: erf(2) erf(2) sage: erf(2).n() 0.995322265018953 sage: erf(2).n(100) 0.99532226501895273416206925637 sage: erf(ComplexField(100)(2+3j)) -20.829461427614568389103088452 + 8.6873182714701631444280787545*I Basic symbolic properties are handled by Sage and Maxima:: sage: x = var("x") sage: diff(erf(x),x) 2*e^(-x^2)/sqrt(pi) sage: integrate(erf(x),x) x*erf(x) + e^(-x^2)/sqrt(pi) ALGORITHM: Numerical evaluation is handled using mpmath, but symbolics are handled by Sage and Maxima. REFERENCES: - http://en.wikipedia.org/wiki/Error_function - mpmath documentation: error-functions_ TESTS: Check limits:: sage: limit(erf(x),x=0) 0 sage: limit(erf(x),x=infinity) 1 Check that it's odd:: sage: erf(1.0) 0.842700792949715 sage: erf(-1.0) -0.842700792949715 Check against other implementations and against the definition:: sage: erf(3).n() 0.999977909503001 sage: maxima.erf(3).n() 0.999977909503001 sage: (1-pari(3).erfc()) 0.999977909503001 sage: RR(3).erf() 0.999977909503001 sage: (integrate(exp(-x**2),(x,0,3))*2/sqrt(pi)).n() 0.999977909503001 Make sure big numbers don't crash it:: sage: set(erf(45*10**i).n() for i in range(10)) set([1.00000000000000]) Trac #9044:: sage: N(erf(sqrt(2)),200) 0.95449973610364158559943472566693312505644755259664313203267 Trac #11626:: sage: n(erf(2),100) 0.99532226501895273416206925637 sage: erf(2).n(100) 0.99532226501895273416206925637 Test (indirectly) trac #11885:: sage: erf(float(0.5)) 0.52049987781304652 sage: erf(complex(0.5)) (0.52049987781304652+0j) Ensure conversion from maxima elements works:: sage: merf = maxima(erf(x)).sage().operator() sage: merf == erf True Make sure we can dump and load it:: sage: loads(dumps(erf(2))) erf(2) Special-case 0 for immediate evaluation:: sage: erf(0) 0 Make sure that we can hold:: sage: erf(0,hold=True) erf(0) sage: simplify(erf(0,hold=True)) 0 Check that high-precision ComplexField inputs don't trigger an mpmath OverflowError (strange bug while developing):: sage: CC(erf(ComplexField(1000)(2+3j))) -20.8294614276146 + 8.68731827147016*I """ def __init__(self): r""" The error function, defined as \text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt. Sage currently only implements the error function (via a call to PARI) when the input is real. """ See the docstring for :meth:Function_erf. EXAMPLES:: sage: erf(2) erf(2) sage: erf(2).n() sage: erf(1) erf(1) """ BuiltinFunction.__init__(self, "erf", nargs=1, latex_name=r'{\rm erf}', conversions=dict(maxima='erf')) def _eval_(self, x): """ EXAMPLES:: sage: erf(1) erf(1) sage: erf(2.0) 0.995322265018953 sage: loads(dumps(erf)) erf The following fails because we haven't implemented erf yet for complex values:: sage: complex(erf(3*I)) Traceback (most recent call last): ... TypeError: unable to simplify to complex approximation sage: CC(erf(3*I)) 1629.99462260157*I TESTS: Check if conversion from maxima elements work:: Test that trac #8983 is fixed -- although it appears to have been fixed by maxima-related changes, let's keep it that way:: sage: merf = maxima(erf(x)).sage().operator() sage: merf == erf True sage: erf(0) 0 sage: solve(erf(x)==0,x) [x == 0] Verify we're returning the appropriate zero:: sage: erf(0) 0 sage: erf(0.0) 0.000000000000000 sage: erf(RealField(100)(0)) 0.00000000000000000000000000000 """ BuiltinFunction.__init__(self, "erf", latex_name=r"\text{erf}") x_zero = False if not isinstance(x, Expression): if is_inexact(x): return self._evalf_(x, parent(x)) if not x: x_zero = True else: if x._is_numerically_zero(): x_zero = True if x_zero: return parent(x)(0) return None # leaves the expression unevaluated def _evalf_(self, x, parent=None): """ sage: erf(2).n() 0.995322265018953 sage: erf(2).n(150) Traceback (most recent call last): ... NotImplementedError: erf not implemented for precision higher than 53 0.99532226501895273416206925636725292861089180 sage: erf(ComplexField(150)(2j)) 18.564802414575552598704291913241017198858002*I """ try: prec = parent.prec() except AttributeError: # not a Sage parent prec = 0 if prec > 53: raise NotImplementedError, "erf not implemented for precision higher than 53" return parent(1 - pari(float(x)).erfc()) if hasattr(x,'erf'): return x.erf() if isinstance(x, float): return float(RDF(x).erf()) import mpmath return mpmath_utils.call(mpmath.erf, x, parent=parent) def _derivative_(self, x, diff_param=None): """ Derivative of erf function sage: erf(x).diff(x) 2*e^(-x^2)/sqrt(pi) TESTS:: TESTS: Check if #8568 is fixed::