# HG changeset patch
# User Jeroen Demeyer <jdemeyer@cage.ugent.be>
# Date 1319474587 7200
# Node ID 5da8e8e0d10e0abc38bcd88f5099044d5370cdfb
# Parent bda8f14b0bf36072c1a064b2acdf17a94c78d1ee
Fix numeric evaluation of error function
diff git a/sage/functions/other.py b/sage/functions/other.py
a

b


28  28  _eval_ = BuiltinFunction._eval_default 
29  29  def __init__(self): 
30  30  r""" 
31   The error function, defined as 
 31  The error function, defined for real values as 
32  32  `\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} dt`. 
 33  This function is also defined for complex values, via analytic 
 34  continuation. 
33  35  
34   Sage currently only implements the error function (via a call to 
35   PARI) when the input is real. 
 36  Sage implements the error function via the ``erfc()`` function in PARI. 
36  37  
37  38  EXAMPLES:: 
38  39  
… 
… 

43  44  sage: loads(dumps(erf)) 
44  45  erf 
45  46  
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 
53   
54  47  TESTS: 
55  48  
56  49  Check if conversion from maxima elements work:: 
… 
… 

61  54  """ 
62  55  BuiltinFunction.__init__(self, "erf", latex_name=r"\text{erf}") 
63  56  
64   def _evalf_(self, x, parent=None): 
 57  def _evalf_(self, x, parent): 
65  58  """ 
66  59  EXAMPLES:: 
67  60  
68  61  sage: erf(2).n() 
69  62  0.995322265018953 
70   sage: erf(2).n(150) 
71   Traceback (most recent call last): 
72   ... 
73   NotImplementedError: erf not implemented for precision higher than 53 
 63  sage: erf(2).n(200) 
 64  0.99532226501895273416206925636725292861089179704006007673835 
 65  sage: erf(pi  1/2*I).n(100) 
 66  1.0000111669099367825726058952 + 1.6332655417638522934072124548e6*I 
 67  
 68  TESTS: 
 69  
 70  Check that PARI/GP through the GP interface gives the same answer:: 
 71  
 72  sage: gp.set_real_precision(59) # random 
 73  38 
 74  sage: print gp.eval("1  erfc(1)"); print erf(1).n(200); 
 75  0.84270079294971486934122063508260925929606699796630290845994 
 76  0.84270079294971486934122063508260925929606699796630290845994 
74  77  """ 
75  78  try: 
76  79  prec = parent.prec() 
77  80  except AttributeError: # not a Sage parent 
78  81  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  return parent(1)  parent(pari(x).erfc(prec)) 
82  83  
83  84  def _derivative_(self, x, diff_param=None): 
84  85  """ 
diff git a/sage/symbolic/expression.pyx b/sage/symbolic/expression.pyx
a

b


988  988  EXAMPLES:: 
989  989  
990  990  sage: complex(I) 
991   1j 
 991  1j 
992  992  sage: complex(erf(3*I)) 
993   Traceback (most recent call last): 
994   ... 
995   TypeError: unable to simplify to complex approximation 
 993  (1.0000000000000002+1629.8673238578601j) 
996  994  """ 
997  995  try: 
998  996  return self._eval_self(complex) 
… 
… 

2927  2925  sage: taylor(a*log(z), z, 2, 3) 
2928  2926  1/24*(z  2)^3*a  1/8*(z  2)^2*a + 1/2*(z  2)*a + a*log(2) 
2929  2927  
2930   :: 
 2928  :: 
2931  2929  
2932  2930  sage: taylor(sqrt (sin(x) + a*x + 1), x, 0, 3) 
2933  2931  1/48*(3*a^3 + 9*a^2 + 9*a  1)*x^3  1/8*(a^2 + 2*a + 1)*x^2 + 1/2*(a + 1)*x + 1 
2934  2932  
2935   :: 
 2933  :: 
2936  2934  
2937  2935  sage: taylor (sqrt (x + 1), x, 0, 5) 
2938  2936  7/256*x^5  5/128*x^4 + 1/16*x^3  1/8*x^2 + 1/2*x + 1 
2939  2937  
2940   :: 
 2938  :: 
2941  2939  
2942  2940  sage: taylor (1/log (x + 1), x, 0, 3) 
2943  2941  19/720*x^3 + 1/24*x^2  1/12*x + 1/x + 1/2 
2944  2942  
2945   :: 
 2943  :: 
2946  2944  
2947  2945  sage: taylor (cos(x)  sec(x), x, 0, 5) 
2948  2946  1/6*x^4  x^2 
2949  2947  
2950   :: 
 2948  :: 
2951  2949  
2952  2950  sage: taylor ((cos(x)  sec(x))^3, x, 0, 9) 
2953  2951  1/2*x^8  x^6 
2954  2952  
2955   :: 
 2953  :: 
2956  2954  
2957  2955  sage: taylor (1/(cos(x)  sec(x))^3, x, 0, 5) 
2958  2956  15377/7983360*x^4  6767/604800*x^2 + 11/120/x^2 + 1/2/x^4  1/x^6  347/15120 
2959  2957  
2960   
2961   Ticket #7472 fixed (Taylor polynomial in more variables) :: 
 2958  TESTS: 
 2959  
 2960  Check that ticket #7472 is fixed (Taylor polynomial in more variables):: 
2962  2961  
2963  2962  sage: x,y=var('x y'); taylor(x*y^3,(x,1),(y,1),4) 
2964  2963  (y  1)^3*(x  1) + (y  1)^3 + 3*(y  1)^2*(x  1) + 3*(y  1)^2 + 3*(y  1)*(x  1) + x + 3*y  3 
… 
… 

7548  7547   ``multiplicities``  bool (default: False); if True, 
7549  7548  return corresponding multiplicities. This keyword is 
7550  7549  incompatible with ``to_poly_solve=True`` and does not make 
7551   any sense when solving inequality. 
 7550  any sense when solving inequality. 
7552  7551  
7553  7552   ``solution_dict``  bool (default: False); if True or nonzero, 
7554  7553  return a list of dictionaries containing solutions. Not used 
… 
… 

7562  7561  Maxima's ``to_poly_solver`` package to search for more possible 
7563  7562  solutions, but possibly encounter approximate solutions. 
7564  7563  This keyword is incompatible with ``multiplicities=True`` 
7565   and is not used when solving inequality. Setting ``to_poly_solve`` 
 7564  and is not used when solving inequality. Setting ``to_poly_solve`` 
7566  7565  to 'force' (string) omits Maxima's solve command (usefull when 
7567  7566  some solution of trigonometric equations are lost). 
7568  7567  
… 
… 

7610  7609  sage: solve(Q*sqrt(Q^2 + 2)  1, Q) 
7611  7610  [Q == 1/sqrt(Q^2 + 2)] 
7612  7611  sage: solve(Q*sqrt(Q^2 + 2)  1, Q, to_poly_solve=True) 
7613   [Q == 1/sqrt(sqrt(2) + 1), Q == 1/sqrt(sqrt(2) + 1)] 
 7612  [Q == 1/sqrt(sqrt(2) + 1), Q == 1/sqrt(sqrt(2) + 1)] 
7614  7613  
7615  7614  In some cases there may be infinitely many solutions indexed 
7616  7615  by a dummy variable. If it begins with ``z``, it is implicitly 