Ticket #13050: trac_13050_erf_mpmath_eval.patch

File trac_13050_erf_mpmath_eval.patch, 5.1 KB (added by benjaminfjones, 8 years ago)

add mpmath evaluation to erf

  • sage/functions/other.py

    # HG changeset patch
    # User Benjamin Jones <benjaminfjones@gmail.com>
    # Date 1338226746 25200
    # Node ID 1ec071c0fbe203a01691db0c827fd5c0a928573e
    # Parent  c696547b633b9e159ce921b1aab575b57a803cb0
    Trac 13050: add mpmath evaluation to erf
    
    diff --git a/sage/functions/other.py b/sage/functions/other.py
    a b  
    1414
    1515import sage.structure.element
    1616coercion_model = sage.structure.element.get_coercion_model()
     17from sage.structure.coerce import parent as sage_structure_coerce_parent
    1718
    18 from sage.structure.coerce import parent
    1919from sage.symbolic.constants import pi
    2020from sage.symbolic.function import is_inexact
    2121from sage.functions.log import exp
     
    4646            sage: erf(2).n(200)
    4747            0.99532226501895273416206925636725292861089179704006007673835
    4848
     49        The function supports numerical evaluation using mpmath (default) or PARI::
     50
     51            sage: erf(2).n(algorithm='pari')
     52            0.995322265018953
     53            sage: erf(2).n(200, algorithm='pari')
     54            0.99532226501895273416206925636725292861089179704006007673835
     55
     56        WARNING: PARI may overflow if the input to `erf` is too large::
     57
     58            sage: erf(10**10).n(algorithm='pari')
     59            ...
     60            PariError:  (15)
     61
     62        The mpmath library handles large numbers gracefully::
     63
     64            sage: erf(10**10).n(algorithm='mpmath')
     65            1.00000000000000
     66
    4967        Complex values can be evaluated as well::
    5068
    5169            sage: erf(3*I).n()
     
    101119        """
    102120        if not isinstance(x, Expression):
    103121            if is_inexact(x):
    104                 return self._evalf_(x, parent=parent(x))
     122                return self._evalf_(x, parent=sage_structure_coerce_parent(x))
    105123            elif x == Integer(0):
    106124                return Integer(0)
    107125        elif x.is_trivial_zero():
     
    114132
    115133            sage: erf(2).n()
    116134            0.995322265018953
     135            sage: erf(2).n(algorithm='pari')
     136            0.995322265018953
    117137            sage: erf(2).n(200)
    118138            0.99532226501895273416206925636725292861089179704006007673835
    119139            sage: erf(pi - 1/2*I).n(100)
     140            1.0000111669099367825726058952 + 1.6332655417638522934072124547e-6*I
     141            sage: erf(pi - 1/2*I).n(100, algorithm='pari')
    120142            1.0000111669099367825726058952 + 1.6332655417638522934072124548e-6*I
    121143
    122144        TESTS:
     
    125147
    126148            sage: gp.set_real_precision(59)  # random
    127149            38
    128             sage: print gp.eval("1 - erfc(1)"); print erf(1).n(200);
     150            sage: print gp.eval("1 - erfc(1)"); print erf(1).n(200, algorithm='pari');
    129151            0.84270079294971486934122063508260925929606699796630290845994
    130152            0.84270079294971486934122063508260925929606699796630290845994
    131153        """
    132         try:
    133             prec = parent.prec()
    134         except AttributeError: # not a Sage parent
    135             prec = 0
    136         return parent(1) - parent(pari(x).erfc(prec))
    137        
     154        if not algorithm or algorithm == 'mpmath': # default algorithm is 'mpmath'
     155            R = parent or sage_structure_coerce_parent(x)
     156            import mpmath
     157            return mpmath_utils.call(mpmath.erf, x, parent=R)
     158        elif algorithm == 'pari':
     159            try:
     160                prec = parent.prec()
     161            except AttributeError: # not a Sage parent
     162                prec = 0
     163            return parent(1) - parent(pari(x).erfc(prec))
     164        else:
     165            raise NotImplementedError("erf only supports algorithm = 'mpmath' (default), "\
     166                                      "or 'pari'")
     167
    138168    def _derivative_(self, x, diff_param=None):
    139169        """
    140170        Derivative of erf function
     
    850878        if not isinstance(x, Expression) and not isinstance(y, Expression) and \
    851879               (is_inexact(x) or is_inexact(y)):
    852880            x, y = coercion_model.canonical_coercion(x, y)
    853             return self._evalf_(x, y, parent(x))
     881            return self._evalf_(x, y, sage_structure_coerce_parent(x))
    854882
    855883        if y == 0:
    856884            return gamma(x)
     
    16711699
    16721700        """
    16731701        if not isinstance(x,Expression): # x contains no variables
    1674             if parent(x)(0)==x: #compatibility with maxima
    1675                 return parent(x)(0)
     1702            if sage_structure_coerce_parent(x)(0)==x: #compatibility with maxima
     1703                return sage_structure_coerce_parent(x)(0)
    16761704            else:
    16771705                if is_inexact(x): # inexact complex numbers, e.g. 2.0+i
    1678                     return self._evalf_(x, parent(x))
     1706                    return self._evalf_(x, sage_structure_coerce_parent(x))
    16791707                else:  # exact complex numbers, e.g. 2+i
    16801708                    return arctan2(imag_part(x),real_part(x))
    16811709        else:
     
    17161744            pass
    17171745        # try to find a parent that support .arg()
    17181746        if parent_d is None:
    1719             parent_d = parent(x)
     1747            parent_d = sage_structure_coerce_parent(x)
    17201748        try:
    17211749            parent_d = parent_d.complex_field()
    17221750        except AttributeError: