Ticket #1173: trac_1173_add_complex_argument_erf.patch

File trac_1173_add_complex_argument_erf.patch, 5.6 KB (added by dsm, 6 years ago)
  • sage/functions/other.py

    # HG changeset patch
    # User D. S. McNeil
    # Date 1298514253 -28800
    # Node ID 678841360c1b8902336e2138ca893236bbc5d7cf
    # Parent  8438b7c20d79c02a2ece3e1c3f7224a772ff8f07
    Trac 1173: add complex-argument erf
    
    diff -r 8438b7c20d79 -r 678841360c1b sage/functions/other.py
    a b  
    2424one_half = ~SR(2)
    2525
    2626class Function_erf(BuiltinFunction):
    27     _eval_ = BuiltinFunction._eval_default
    2827    def __init__(self):
    2928        r"""
    3029        The error function, defined as
    3130        `\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt`.
    3231
    33         Sage currently only implements the error function (via a call to
    34         PARI) when the input is real.
     32        Sage implements the error function for both real and complex arguments
     33        (via mpmath).
    3534
    3635        EXAMPLES::
    37 
     36       
    3837            sage: erf(2)
    3938            erf(2)
    4039            sage: erf(2).n()
    4140            0.995322265018953
    42             sage: loads(dumps(erf))
    43             erf
     41            sage: erf(2).n(100)
     42            0.99532226501895273416206925637
     43            sage: erf(ComplexField(100)(2+3j))
     44            -20.829461427614568389103088452 + 8.6873182714701631444280787545*I
     45            sage: a = sqrt(pi)*I*erf(2*I)/2
     46            sage: a
     47            1/2*I*sqrt(pi)*erf(2*I)
     48            sage: CC(a)
     49            -16.4526277655072
     50           
     51        Note that erf(0) is immediately evaluated, but this behaviour
     52        can be suppressed via the hold parameter, and then undone
     53        via simplify:
    4454
    45         The following fails because we haven't implemented
    46         erf yet for complex values::
    47        
    48             sage: complex(erf(3*I))
    49             Traceback (most recent call last):
    50             ...
    51             TypeError: unable to simplify to complex approximation
     55            sage: erf(0)
     56            0
     57            sage: erf(0,hold=True)
     58            erf(0)
     59            sage: simplify(erf(0,hold=True))
     60            0
    5261
    53         TESTS:
     62
     63        TESTS::
    5464
    5565        Check if conversion from maxima elements work::
    5666
    5767            sage: merf = maxima(erf(x)).sage().operator()
    5868            sage: merf == erf
    5969            True
     70
     71        Make sure we can dump and load it:
     72
     73            sage: loads(dumps(erf(2)))
     74            erf(2)
     75
     76        Check to make sure that exact ring inputs return unevaluated
     77        (note that zero is an exception), that inexact rings evaluate,
     78        and that everything returns an appropriate type:
     79
     80            sage: exact = SR, ZZ, QQ, int, long
     81            sage: list(erf(ring(1)) for ring in exact)
     82            [erf(1), erf(1), erf(1), erf(1), erf(1L)]
     83            sage: erf(e), erf(pi)
     84            (erf(e), erf(pi))
     85            sage: inexact = (RR, RealField(100), RDF, float, CC, ComplexField(100), CDF, complex)
     86            sage: list(erf(ring(1)) for ring in inexact)
     87            [0.842700792949715, 0.84270079294971486934122063508, 0.84270079295, 0.84270079294971489,
     88            0.842700792949715, 0.84270079294971486934122063508, 0.84270079295, (0.84270079294971489+0j)]
     89            sage: list(parent(ring(1)) == parent(erf(ring(1))) for ring in inexact)
     90            [True, True, True, True, True, True, True, True]
     91
     92        Make sure that we can hold:
     93
     94            sage: erf(ComplexField(100)(2+3j),hold=True)
     95            erf(2.0000000000000000000000000000 + 3.0000000000000000000000000000*I)
     96           
     97        Check that large-precision ComplexField inputs don't trigger an mpmath OverflowError
     98        (strange bug while developing):
     99
     100            sage: CC(erf(ComplexField(1000)(2+3j)))
     101            -20.8294614276146 + 8.68731827147016*I
     102           
    60103        """
    61104        BuiltinFunction.__init__(self, "erf", latex_name=r"\text{erf}")
    62105
     106    def _eval_(self, x):
     107        """
     108        EXAMPLES::
     109
     110            sage: erf(1)
     111            erf(1)
     112            sage: erf(2.0)
     113            0.995322265018953
     114
     115        TESTS::
     116
     117        Test that trac #8983 is fixed -- although it appears to have been fixed
     118        by maxima-related changes, let's keep it that way:
     119
     120            sage: erf(0)
     121            0
     122            sage: solve(erf(x)==0,x)
     123            [x == 0]
     124
     125        Verify we're returning the appropriate zero:
     126
     127            sage: erf(0)
     128            0
     129            sage: erf(0.0)
     130            0.000000000000000
     131            sage: erf(RealField(100)(0))
     132            0.00000000000000000000000000000
     133       
     134        """
     135        if not isinstance(x, Expression) and is_inexact(x):
     136            return self._evalf_(x, parent(x))
     137        if x == 0:
     138            return parent(x)(0)
     139
    63140    def _evalf_(self, x, parent=None):
    64141        """
    65142        EXAMPLES::
     
    67144            sage: erf(2).n()
    68145            0.995322265018953
    69146            sage: erf(2).n(150)
    70             Traceback (most recent call last):
    71             ...
    72             NotImplementedError: erf not implemented for precision higher than 53
     147            0.99532226501895273416206925636725292861089180
     148            sage: erf(ComplexField(150)(2j))
     149            18.564802414575552598704291913241017198858002*I
     150
    73151        """
     152        from mpmath import erf
    74153        try:
    75154            prec = parent.prec()
    76         except AttributeError: # not a Sage parent
    77             prec = 0
    78         if prec > 53:
    79             raise NotImplementedError, "erf not implemented for precision higher than 53"
    80         return parent(1 - pari(float(x)).erfc())
    81        
     155        except AttributeError:
     156            prec = 53
     157        return parent(mpmath_utils.call(erf, x, prec=prec))
     158 
    82159    def _derivative_(self, x, diff_param=None):
    83160        """
    84161        Derivative of erf function