Ticket #1173: trac_1173_complex_erf_v3.patch

File trac_1173_complex_erf_v3.patch, 7.4 KB (added by dsm, 8 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 b  
    77from sage.symbolic.pynac import py_factorial_py
    88from sage.libs.pari.gen import pari
    99from sage.symbolic.all import SR
    10 from sage.rings.all import Integer, Rational, RealField, RR, \
     10from sage.rings.all import Integer, Rational, RealField, RR, RDF, \
    1111     is_ComplexNumber, ComplexField
    1212from sage.misc.latex import latex
    1313import math
    … …  
    2424
    2525one_half = ~SR(2)
    2626
     27
    2728class 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: `error-functions`_
     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: (1-pari(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    Special-case 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 high-precision 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    """
    29153    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`.
    36156
    37157        EXAMPLES::
    38158
    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)
    42173            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
    53176
    54177        TESTS:
    55178
    56         Check if conversion from maxima elements work::
     179        Test that trac #8983 is fixed -- although it appears to have been fixed
     180        by maxima-related changes, let's keep it that way::
    57181
    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
    61197        """
    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
    63211
    64212    def _evalf_(self, x, parent=None):
    65213        """
    … …  
    68216            sage: erf(2).n()
    69217            0.995322265018953
    70218            sage: erf(2).n(150)
    71             Traceback (most recent call last):
    72             ...
    73             NotImplementedError: erf not implemented for precision higher than 53
     219            0.99532226501895273416206925636725292861089180
     220            sage: erf(ComplexField(150)(2j))
     221            18.564802414575552598704291913241017198858002*I
    74222        """
    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
    83230    def _derivative_(self, x, diff_param=None):
    84231        """
    85232        Derivative of erf function
    … …  
    89236            sage: erf(x).diff(x)
    90237            2*e^(-x^2)/sqrt(pi)
    91238
    92         TESTS::
     239        TESTS:
    93240
    94241        Check if #8568 is fixed::
    95242