Ticket #1173: trac_1173_complex_erf_v3.patch
File trac_1173_complex_erf_v3.patch, 7.4 KB (added by , 7 years ago) 


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 b 7 7 from sage.symbolic.pynac import py_factorial_py 8 8 from sage.libs.pari.gen import pari 9 9 from sage.symbolic.all import SR 10 from sage.rings.all import Integer, Rational, RealField, RR, \10 from sage.rings.all import Integer, Rational, RealField, RR, RDF, \ 11 11 is_ComplexNumber, ComplexField 12 12 from sage.misc.latex import latex 13 13 import math … … 24 24 25 25 one_half = ~SR(2) 26 26 27 27 28 class Function_erf(BuiltinFunction): 28 _eval_ = BuiltinFunction._eval_default 29 r""" 30 The error function, defined for real values as 31 32 `\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} dt`. 33 34 This function is also defined for complex values, via analytic 35 continuation. 36 37 38 EXAMPLES: 39 40 We can evaluate numerically via mpmath:: 41 42 sage: erf(2) 43 erf(2) 44 sage: erf(2).n() 45 0.995322265018953 46 sage: erf(2).n(100) 47 0.99532226501895273416206925637 48 sage: erf(ComplexField(100)(2+3j)) 49 20.829461427614568389103088452 + 8.6873182714701631444280787545*I 50 51 Basic symbolic properties are handled by Sage and Maxima:: 52 53 sage: x = var("x") 54 sage: diff(erf(x),x) 55 2*e^(x^2)/sqrt(pi) 56 sage: integrate(erf(x),x) 57 x*erf(x) + e^(x^2)/sqrt(pi) 58 59 ALGORITHM: 60 61 Numerical evaluation is handled using mpmath, but symbolics are handled 62 by Sage and Maxima. 63 64 REFERENCES: 65 66  http://en.wikipedia.org/wiki/Error_function 67  mpmath documentation: `errorfunctions`_ 68 69 70 TESTS: 71 72 Check limits:: 73 74 sage: limit(erf(x),x=0) 75 0 76 sage: limit(erf(x),x=infinity) 77 1 78 79 Check that it's odd:: 80 81 sage: erf(1.0) 82 0.842700792949715 83 sage: erf(1.0) 84 0.842700792949715 85 86 Check against other implementations and against the definition:: 87 88 sage: erf(3).n() 89 0.999977909503001 90 sage: maxima.erf(3).n() 91 0.999977909503001 92 sage: (1pari(3).erfc()) 93 0.999977909503001 94 sage: RR(3).erf() 95 0.999977909503001 96 sage: (integrate(exp(x**2),(x,0,3))*2/sqrt(pi)).n() 97 0.999977909503001 98 99 Make sure big numbers don't crash it:: 100 101 sage: set(erf(45*10**i).n() for i in range(10)) 102 set([1.00000000000000]) 103 104 Trac #9044:: 105 106 sage: N(erf(sqrt(2)),200) 107 0.95449973610364158559943472566693312505644755259664313203267 108 109 Trac #11626:: 110 111 sage: n(erf(2),100) 112 0.99532226501895273416206925637 113 sage: erf(2).n(100) 114 0.99532226501895273416206925637 115 116 Test (indirectly) trac #11885:: 117 118 sage: erf(float(0.5)) 119 0.52049987781304652 120 sage: erf(complex(0.5)) 121 (0.52049987781304652+0j) 122 123 Ensure conversion from maxima elements works:: 124 125 sage: merf = maxima(erf(x)).sage().operator() 126 sage: merf == erf 127 True 128 129 Make sure we can dump and load it:: 130 131 sage: loads(dumps(erf(2))) 132 erf(2) 133 134 Specialcase 0 for immediate evaluation:: 135 136 sage: erf(0) 137 0 138 139 Make sure that we can hold:: 140 141 sage: erf(0,hold=True) 142 erf(0) 143 sage: simplify(erf(0,hold=True)) 144 0 145 146 Check that highprecision ComplexField inputs don't trigger 147 an mpmath OverflowError (strange bug while developing):: 148 149 sage: CC(erf(ComplexField(1000)(2+3j))) 150 20.8294614276146 + 8.68731827147016*I 151 152 """ 29 153 def __init__(self): 30 r""" 31 The error function, defined as 32 `\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} dt`. 33 34 Sage currently only implements the error function (via a call to 35 PARI) when the input is real. 154 """ 155 See the docstring for :meth:`Function_erf`. 36 156 37 157 EXAMPLES:: 38 158 39 sage: erf(2) 40 erf(2) 41 sage: erf(2).n() 159 sage: erf(1) 160 erf(1) 161 162 """ 163 BuiltinFunction.__init__(self, "erf", nargs=1, latex_name=r'{\rm erf}', 164 conversions=dict(maxima='erf')) 165 166 def _eval_(self, x): 167 """ 168 EXAMPLES:: 169 170 sage: erf(1) 171 erf(1) 172 sage: erf(2.0) 42 173 0.995322265018953 43 sage: loads(dumps(erf)) 44 erf 45 46 The following fails because we haven't implemented 47 erf yet for complex values:: 48 49 sage: complex(erf(3*I)) 50 Traceback (most recent call last): 51 ... 52 TypeError: unable to simplify to complex approximation 174 sage: CC(erf(3*I)) 175 1629.99462260157*I 53 176 54 177 TESTS: 55 178 56 Check if conversion from maxima elements work:: 179 Test that trac #8983 is fixed  although it appears to have been fixed 180 by maximarelated changes, let's keep it that way:: 57 181 58 sage: merf = maxima(erf(x)).sage().operator() 59 sage: merf == erf 60 True 182 sage: erf(0) 183 0 184 sage: solve(erf(x)==0,x) 185 [x == 0] 186 187 Verify we're returning the appropriate zero:: 188 189 sage: erf(0) 190 0 191 sage: erf(0.0) 192 0.000000000000000 193 sage: erf(RealField(100)(0)) 194 0.00000000000000000000000000000 195 196 61 197 """ 62 BuiltinFunction.__init__(self, "erf", latex_name=r"\text{erf}") 198 x_zero = False 199 if not isinstance(x, Expression): 200 if is_inexact(x): 201 return self._evalf_(x, parent(x)) 202 if not x: 203 x_zero = True 204 else: 205 if x._is_numerically_zero(): 206 x_zero = True 207 if x_zero: 208 return parent(x)(0) 209 210 return None # leaves the expression unevaluated 63 211 64 212 def _evalf_(self, x, parent=None): 65 213 """ … … 68 216 sage: erf(2).n() 69 217 0.995322265018953 70 218 sage: erf(2).n(150) 71 Traceback (most recent call last):72 ...73 NotImplementedError: erf not implemented for precision higher than 53219 0.99532226501895273416206925636725292861089180 220 sage: erf(ComplexField(150)(2j)) 221 18.564802414575552598704291913241017198858002*I 74 222 """ 75 try: 76 prec = parent.prec() 77 except AttributeError: # not a Sage parent 78 prec = 0 79 if prec > 53: 80 raise NotImplementedError, "erf not implemented for precision higher than 53" 81 return parent(1  pari(float(x)).erfc()) 82 223 if hasattr(x,'erf'): 224 return x.erf() 225 if isinstance(x, float): 226 return float(RDF(x).erf()) 227 import mpmath 228 return mpmath_utils.call(mpmath.erf, x, parent=parent) 229 83 230 def _derivative_(self, x, diff_param=None): 84 231 """ 85 232 Derivative of erf function … … 89 236 sage: erf(x).diff(x) 90 237 2*e^(x^2)/sqrt(pi) 91 238 92 TESTS: :239 TESTS: 93 240 94 241 Check if #8568 is fixed:: 95 242