Repository: sympy/sympy
Pull request: Support Ei(z) expansion at -oo
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #26916
Closes #26937 

#### Brief description of what is fixed or changed
The series expansion of `Ei(x)` at `oo` and `-oo` is the same hence we just had to tweak the condition

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* series
  * Support Ei(z) expansion at -oo
<!-- END RELEASE NOTES -->


Problem statement:
Ei(-x) * exp(x) (exponential integral) limit to infinity error out
The following limit computes fine:

```
(Ei(x) * exp(-x)).limit(x, oo)
# Result: 0
```

But the following doesn't:

```
(Ei(-x) * exp(x)).limit(x, oo)
# Or equivalent: (Ei(x) * exp(-x)).limit(x, -oo)
# Expected: 0
```

Instead, sympy will error out with

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~\AppData\Roaming\Python\Python39\site-packages\sympy\series\gruntz.py:557, in mrv_leadterm(e, x)
    556 try:
--> 557     lt = f.leadterm(w, logx=logw)
    558 except (NotImplementedError, PoleError, ValueError):

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\core\expr.py:3533, in Expr.leadterm(self, x, logx, cdir)
   3532 if x in c.free_symbols:
-> 3533     raise ValueError(filldedent("""
   3534         cannot compute leadterm(%s, %s). The coefficient
   3535         should have been free of %s but got %s""" % (self, x, x, c)))
   3536 c = c.subs(d, log(x))

ValueError: 
cannot compute leadterm(_eis(-1/_w), _w). The coefficient should have
been free of _w but got _eis(-1/_w)

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[84], line 2
      1 (Ei(x) * exp(-x)).limit(x, oo)
----> 2 (Ei(x) * exp(-x)).limit(x, -oo)

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\core\expr.py:3418, in Expr.limit(self, x, xlim, dir)
   3415 """ Compute limit x->xlim.
   3416 """
   3417 from sympy.series.limits import limit
-> 3418 return limit(self, x, xlim, dir)

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\series\limits.py:64, in limit(e, z, z0, dir)
     13 def limit(e, z, z0, dir="+"):
     14     """Computes the limit of ``e(z)`` at the point ``z0``.
     15 
     16     Parameters
   (...)
     61      limit_seq : returns the limit of a sequence.
     62     """
---> 64     return Limit(e, z, z0, dir).doit(deep=False)

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\series\limits.py:375, in Limit.doit(self, **hints)
    372 l = None
    374 try:
--> 375     r = gruntz(e, z, z0, dir)
    376     if r is S.NaN or l is S.NaN:
    377         raise PoleError()

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\series\gruntz.py:733, in gruntz(e, z, z0, dir)
    730     else:
    731         raise NotImplementedError("dir must be '+' or '-'")
--> 733 r = limitinf(e0, z)
    735 # This is a bit of a heuristic for nice results... we always rewrite
    736 # tractable functions in terms of familiar intractable ones.
    737 # It might be nicer to rewrite the exactly to what they were initially,
    738 # but that would take some work to implement.
    739 return r.rewrite('intractable', deep=True)

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\core\cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
     69 @wraps(func)
     70 def wrapper(*args, **kwargs):
     71     try:
---> 72         retval = cfunc(*args, **kwargs)
     73     except TypeError as e:
     74         if not e.args or not e.args[0].startswith('unhashable type:'):

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\series\gruntz.py:453, in limitinf(e, x)
    451     c0, e0 = mrv_leadterm(e.min, x)
    452 else:
--> 453     c0, e0 = mrv_leadterm(e, x)
    454 sig = sign(e0, x)
    455 if sig == 1:

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\core\cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
     69 @wraps(func)
     70 def wrapper(*args, **kwargs):
     71     try:
---> 72         retval = cfunc(*args, **kwargs)
     73     except TypeError as e:
     74         if not e.args or not e.args[0].startswith('unhashable type:'):

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\series\gruntz.py:563, in mrv_leadterm(e, x)
    561 incr = S.One
    562 while _series.is_Order:
--> 563     _series = f._eval_nseries(w, n=n0+incr, logx=logw)
    564     incr *= 2
    565 series = _series.expand().removeO()

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\functions\special\error_functions.py:2782, in _eis._eval_nseries(self, x, n, logx, cdir)
   2780     f = self._eval_rewrite_as_intractable(*self.args)
   2781     return f._eval_nseries(x, n, logx)
-> 2782 return super()._eval_nseries(x, n, logx)

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\core\function.py:690, in Function._eval_nseries(self, x, n, logx, cdir)
    688 a0 = [t.limit(x, 0) for t in a]
    689 if any(t.has(oo, -oo, zoo, nan) for t in a0):
--> 690     return self._eval_aseries(n, args0, x, logx)
    691 # Careful: the argument goes to oo, but only logarithmically so. We
    692 # are supposed to do a power series expansion "around the
    693 # logarithmic term". e.g.
    694 #      f(1+x+log(x))
    695 #     -> f(1+logx) + x*f'(1+logx) + O(x**2)
    696 # where 'logx' is given in the argument
    697 a = [t._eval_nseries(x, n, logx) for t in args]

File ~\AppData\Roaming\Python\Python39\site-packages\sympy\functions\special\error_functions.py:2751, in _eis._eval_aseries(self, n, args0, x, logx)
   2749 from sympy.series.order import Order
   2750 if args0[0] != S.Infinity:
-> 2751     return super(_erfs, self)._eval_aseries(n, args0, x, logx)
   2753 z = self.args[0]
   2754 l = [factorial(k) * (1/z)**(k + 1) for k in range(n)]

TypeError: super(type, obj): obj must be an instance or subtype of type
```

Hints from issue discussion:
The immediate error message results from a bug that is easily fixed:
```diff
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 09279588b6..ba29c4d7d4 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -2748,7 +2748,7 @@ class _eis(Function):
     def _eval_aseries(self, n, args0, x, logx):
         from sympy.series.order import Order
         if args0[0] != S.Infinity:
-            return super(_erfs, self)._eval_aseries(n, args0, x, logx)
+            return super()._eval_aseries(n, args0, x, logx)
 
         z = self.args[0]
         l = [factorial(k) * (1/z)**(k + 1) for k in range(n)]
```
In this case the limit does not compute but does not exit with the error:
```python
In [1]: (Ei(-x) * exp(x)).limit(x, oo)
Out[1]: 
    ⎛ x       ⎞
lim ⎝ℯ ⋅Ei(-x)⎠
x─→∞  
```
Great thanks! A little bit unfortunate that it can't compute the limit however. Luckily WolframAlpha didn't seem to have any issues with it.
The basic bug is easy to fix.

Evaluating the limit would require having an asymptotic expansion for `Ei(x)` at `-oo`:
```python
In [5]: Ei(x).series(x, oo)
Out[5]: 
⎛120   24   6    2    1        ⎛1        ⎞⎞  x
⎜─── + ── + ── + ── + ─ + 1 + O⎜──; x → ∞⎟⎟⋅ℯ 
⎜ 5     4    3    2   x        ⎜ 6       ⎟⎟   
⎝x     x    x    x             ⎝x        ⎠⎠   
──────────────────────────────────────────────
                      x                       

In [6]: Ei(x).series(x, -oo)
...
PoleError: 
Asymptotic expansion of Ei around [-oo] is not implemented.
```
There is an asymptotic description on [Wikipedia about Exponential Integrals](https://en.wikipedia.org/wiki/Exponential_integral) which is asymptotic beyond all borders in the complex plain. Would that help?
It is possibly just a case of someone needing to implement the formula somewhere.

Commit message: Support Ei(z) expansion at -oo

Patch to review:
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index fd616c49af6f..945bad4c4384 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -1263,7 +1263,7 @@ def _eval_aseries(self, n, args0, x, logx):
         from sympy.series.order import Order
         point = args0[0]
 
-        if point is S.Infinity:
+        if point in (S.Infinity, S.NegativeInfinity):
             z = self.args[0]
             s = [factorial(k) / (z)**k for k in range(n)] + \
                     [Order(1/z**n, x)]
@@ -2766,8 +2766,8 @@ class _eis(Function):
 
     def _eval_aseries(self, n, args0, x, logx):
         from sympy.series.order import Order
-        if args0[0] != S.Infinity:
-            return super(_erfs, self)._eval_aseries(n, args0, x, logx)
+        if args0[0] not in (S.Infinity, S.NegativeInfinity):
+            return super()._eval_aseries(n, args0, x, logx)
 
         z = self.args[0]
         l = [factorial(k) * (1/z)**(k + 1) for k in range(n)]
diff --git a/sympy/functions/special/tests/test_error_functions.py b/sympy/functions/special/tests/test_error_functions.py
index 051e1af4200b..dc699c330b8c 100644
--- a/sympy/functions/special/tests/test_error_functions.py
+++ b/sympy/functions/special/tests/test_error_functions.py
@@ -437,6 +437,11 @@ def test_ei():
     assert Ei(x).series(x, 1, 3) == Ei(1) + E*(x - 1) + O((x - 1)**3, (x, 1))
     assert Ei(x).series(x, oo) == \
         (120/x**5 + 24/x**4 + 6/x**3 + 2/x**2 + 1/x + 1 + O(x**(-6), (x, oo)))*exp(x)/x
+    assert Ei(x).series(x, -oo) == \
+        (120/x**5 + 24/x**4 + 6/x**3 + 2/x**2 + 1/x + 1 + O(x**(-6), (x, -oo)))*exp(x)/x
+    assert Ei(-x).series(x, oo).expand() == \
+        120*exp(-x)/x**6 - 24*exp(-x)/x**5 + 6*exp(-x)/x**4 - 2*exp(-x)/x**3 + \
+        exp(-x)/x**2 - exp(-x)/x + O(exp(-x)/x**7, (x, oo))
 
     assert str(Ei(cos(2)).evalf(n=10)) == '-0.6760647401'
     raises(ArgumentIndexError, lambda: Ei(x).fdiff(2))
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 21777c15e65c..ac28ad3eafbf 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -1412,3 +1412,8 @@ def test_issue_26250():
     e1 = ((1-3*x**2)*e**2/2 - (x**2-2*x+1)*e*k/2)
     e2 = pi**2*(x**8 - 2*x**7 - x**6 + 4*x**5 - x**4 - 2*x**3 + x**2)
     assert limit(e1/e2, x, 0) == -S(1)/8
+
+
+def test_issue_26916():
+    assert limit(Ei(x)*exp(-x), x, +oo) == 0
+    assert limit(Ei(x)*exp(-x), x, -oo) == 0


Code review comments:
diff_hunk: @@ -437,6 +437,11 @@ def test_ei():
     assert Ei(x).series(x, 1, 3) == Ei(1) + E*(x - 1) + O((x - 1)**3, (x, 1))
     assert Ei(x).series(x, oo) == \
         (120/x**5 + 24/x**4 + 6/x**3 + 2/x**2 + 1/x + 1 + O(x**(-6), (x, oo)))*exp(x)/x
+    assert Ei(x).series(x, -oo) == \
+        (120/x**5 + 24/x**4 + 6/x**3 + 2/x**2 + 1/x + 1 + O(x**(-6), (x, -oo)))*exp(x)/x
+    assert Ei(-x).series(x, oo).expand() == \
+        120*exp(-x)/x**6 - 24*exp(-x)/x**5 + 6*exp(-x)/x**4 - 2*exp(-x)/x**3 + \
+        exp(-x)/x**2 - exp(-x)/x + O(exp(-x)/x**7, (x, oo))
path: sympy/functions/special/tests/test_error_functions.py
text: @user1:
Why is `expand` being called in one case but not the other?
```python
In [1]: Ei(-x).series(x, oo)
Out[1]: 
 ⎛  120   24   6    2    1        ⎛1        ⎞⎞  -x 
-⎜- ─── + ── - ── + ── - ─ + 1 + O⎜──; x → ∞⎟⎟⋅ℯ   
 ⎜   5     4    3    2   x        ⎜ 6       ⎟⎟     
 ⎝  x     x    x    x             ⎝x        ⎠⎠     
───────────────────────────────────────────────────
                         x                         

In [2]: Ei(-x).series(x, -oo)
Out[2]: 
 ⎛  120   24   6    2    1        ⎛1         ⎞⎞  -x 
-⎜- ─── + ── - ── + ── - ─ + 1 + O⎜──; x → -∞⎟⎟⋅ℯ   
 ⎜   5     4    3    2   x        ⎜ 6        ⎟⎟     
 ⎝  x     x    x    x             ⎝x         ⎠⎠     
────────────────────────────────────────────────────
                         x                          

In [3]: Ei(-x).series(x, -oo).expand()
Out[3]: 
     -x       -x      -x      -x    -x    -x    ⎛  -x         ⎞
120⋅ℯ     24⋅ℯ     6⋅ℯ     2⋅ℯ     ℯ     ℯ      ⎜-ℯ           ⎟
─────── - ────── + ───── - ───── + ─── - ─── + O⎜─────; x → -∞⎟
   6         5       4       3      2     x     ⎜  7          ⎟
  x         x       x       x      x            ⎝ x           ⎠
```

@author:
Well cause I wasn't able to figure out why something like fails
```
>>> Ei(-x).series(x, oo)
-(-120/x**5 + 24/x**4 - 6/x**3 + 2/x**2 - 1/x + 1 + O(x**(-6), (x, oo)))*exp(-x)/x
>>> -(-120/x**5 + 24/x**4 - 6/x**3 + 2/x**2 - 1/x + 1 + O(x**(-6), (x, oo)))*exp(-x)/x
(120/x**5 - 24/x**4 + 6/x**3 - 2/x**2 + 1/x - 1 + O(x**(-6), (x, oo)))*exp(-x)/x
>>> assert Ei(-x).series(x, oo) == (120/x**5 - 24/x**4 + 6/x**3 - 2/x**2 + 1/x - 1 + O(x**(-6), (x, oo)))*exp(-x)/x
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AssertionError
```

And hence wasn't sure how could I add a test for this and hence I ended up using expand !

The result is as expected. Am I missing something here ?

@user1:
Brackets are needed:
```python
In [3]: -(-120/x**5 + 24/x**4 - 6/x**3 + 2/x**2 - 1/x + 1 + O(x**(-6), (x, oo)))*exp(-x)/x
Out[3]: 
⎛120   24   6    2    1        ⎛1        ⎞⎞  -x
⎜─── - ── + ── - ── + ─ - 1 + O⎜──; x → ∞⎟⎟⋅ℯ  
⎜ 5     4    3    2   x        ⎜ 6       ⎟⎟    
⎝x     x    x    x             ⎝x        ⎠⎠    
───────────────────────────────────────────────
                       x                       

In [4]: -((-120/x**5 + 24/x**4 - 6/x**3 + 2/x**2 - 1/x + 1 + O(x**(-6), (x, oo)))*exp(-x)/x)
Out[4]: 
 ⎛  120   24   6    2    1        ⎛1        ⎞⎞  -x 
-⎜- ─── + ── - ── + ── - ─ + 1 + O⎜──; x → ∞⎟⎟⋅ℯ   
 ⎜   5     4    3    2   x        ⎜ 6       ⎟⎟     
 ⎝  x     x    x    x             ⎝x        ⎠⎠     
───────────────────────────────────────────────────
                         x 
```
This is because of automatic distribution in a 2-arg Mul (a longstanding issue):
```python
In [5]: 2*(x + y)*z
Out[5]: z⋅(2⋅x + 2⋅y)

In [6]: 2*((x + y)*z)
Out[6]: 2⋅z⋅(x + y)
```

@author:
Thanks for the help !

@user1:
Probably in this case distribution is the right thing though and the series code should be changed to return the distributed result.

@author:
Arghh I just committed the change with the brackets :|

@user1:
That's fine. This PR is good. Just a future PR could maybe handle this. Ideally we want the series to be in a well-defined canonical form so that tests like this are not so fragile.