Repository: sympy/sympy
Pull request: Implemented _eval_nseries() for Hyper
<!-- 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 #18193 

#### Brief description of what is fixed or changed
In reference to Issue #18193 the function was implemented. Example of
the working code is as follows:

In [1]: from sympy import hyper

In [2]: from sympy.abc import x, y

In [3]: hyper((1, 2),(1, 2), x)._eval_nseries(x, 7, y)
Out[3]: 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + x**6/720 + O(x**7)

In [4]: hyper((1, 2),(3, 4, 5), x)._eval_nseries(x, 7, y)
Out[4]: 1 + x/30 + x**2/1200 + x**3/63000 + x**4/4233600 + x**5/355622400 + x**6/36578304000 + O(x**7)


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. 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 -->
* functions
   * Added _eval_nseries() functionality to hyper
<!-- END RELEASE NOTES -->

Problem statement:
AttributeError: integrate(1 / (1 + x**4)**(S(1)/4), [x, 0, oo])
```
>>> from sympy import *
>>> x = Symbol('x')
>>> integrate(1 / (1 + x**4)**(S(1)/4), [x, 0, oo])
Traceback (most recent call last):
  File "test.py", line 40, in <module>
    print(integrate(1 / (1 + x**4)**(S(1)/4), [x, 0, oo]))
  File "/home/eward/se/sympy/integrals/integrals.py", line 1522, in integrate
    return integral.doit(**doit_flags)
  File "/home/eward/se/sympy/integrals/integrals.py", line 669, in doit
    evalued = Add(*others)._eval_interval(x, a, b)
  File "/home/eward/se/sympy/core/expr.py", line 964, in _eval_interval
    B = _eval_endpoint(left=False)
  File "/home/eward/se/sympy/core/expr.py", line 949, in _eval_endpoint
    C = limit(self, x, c, "+" if left else "-")
  File "/home/eward/se/sympy/series/limits.py", line 71, in limit
    return Limit(e, z, z0, dir).doit(deep=False)
  File "/home/eward/se/sympy/series/limits.py", line 261, in doit
    r = gruntz(e, z, z0, dir)
  File "/home/eward/se/sympy/series/gruntz.py", line 671, in gruntz
    r = limitinf(e0, z)
  File "/home/eward/se/sympy/core/cache.py", line 94, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/eward/se/sympy/series/gruntz.py", line 435, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "/home/eward/se/sympy/core/cache.py", line 94, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/eward/se/sympy/series/gruntz.py", line 524, in mrv_leadterm
    series = calculate_series(f, w, logx=logw)
  File "/home/eward/se/sympy/series/gruntz.py", line 477, in calculate_series
    for t in e.lseries(x, logx=logx):
  File "/home/eward/se/sympy/core/expr.py", line 2924, in yield_lseries
    for si in s:
  File "/home/eward/se/sympy/core/expr.py", line 3136, in _eval_lseries
    series = self._eval_nseries(x, n=n, logx=logx)
  File "/home/eward/se/sympy/core/mul.py", line 1749, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "/home/eward/se/sympy/core/mul.py", line 1749, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "/home/eward/se/sympy/core/expr.py", line 3223, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "/home/eward/se/sympy/core/function.py", line 673, in _eval_nseries
    a = [t.compute_leading_term(x, logx=logx) for t in args]
  File "/home/eward/se/sympy/core/function.py", line 673, in <listcomp>
    a = [t.compute_leading_term(x, logx=logx) for t in args]
AttributeError: 'TupleArg' object has no attribute 'compute_leading_term'
```

It is trying to compute the `nseries` of a hypergeometric function:

```
>>> hyper((S(1)/4, S(1)/4), (S(5)/4,), exp_polar(I*pi)/x**4).nseries()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/eward/se/sympy/core/expr.py", line 3221, in nseries
    return self.series(x, x0, n, dir)
  File "/home/eward/se/sympy/core/expr.py", line 2870, in series
    rv = self.subs(x, xpos).series(xpos, x0, n, dir, logx=logx)
  File "/home/eward/se/sympy/core/expr.py", line 2877, in series
    s1 = self._eval_nseries(x, n=n, logx=logx)
  File "/home/eward/se/sympy/core/function.py", line 673, in _eval_nseries
    a = [t.compute_leading_term(x, logx=logx) for t in args]
  File "/home/eward/se/sympy/core/function.py", line 673, in <listcomp>
    a = [t.compute_leading_term(x, logx=logx) for t in args]
AttributeError: 'TupleArg' object has no attribute 'compute_leading_term'
```

Hints from issue discussion:
It seems that `hyper` (and other special functions) should have a private `_eval_nseries` method implemented.
I would like to work on this. Could someone please guide where should i start.
I would start by studying the existing implementations. Those can be found by running `git grep 'def _eval_nseries'`.
Why are the following limits taken and how are conditions which decide whether `super()._eval_nseries()` or `self._eval_nseries()` will be called determined ?
https://github.com/sympy/sympy/blob/c31a29a77a8b405d19ef79193ff0878345c62189/sympy/functions/elementary/trigonometric.py#L1136

https://github.com/sympy/sympy/blob/c31a29a77a8b405d19ef79193ff0878345c62189/sympy/functions/special/gamma_functions.py#L195-L196
@jksuom will something like this work:
```
def _eval_nseries(self, x, n, logx, ap, bq):

	from sympy.functions import factorial, RisingFactorial

	x0 = self.args[0].limit(x, 0)

	if not(x0==0):
	return super(gamma, self)._eval_nseries(x, n, logx)

	Rf1 = 1
	Rf2 = 1
	sum = 0

	for i in range(n):
		for p in range(len(ap)):
			Rf1 *= RisingFactorial(a,p+1)
		
		for q in range(len(bq)):
			Rf2 *= RisingFactorial(b,q+1)

		sum += ((Rf1/Rf2)*(x**i))/factorial(i)

	return sum
```
Another idea I had was to use an approach similar to ` _eval_rewrite_as_Sum()` but instead of taking a dummy variable a loop could be implemented
> def _eval_nseries(self, x, n, logx, ap, bq)

`_eval_nseries` has no parameters `ap`, `bq`. If this is for `hyper` (as I assume), then these are the first two arguments of `self`: `self.args[0]` and `self.args[1]`. `x0` will be obtained from `self.args[2]`. This can be seen from the code where instances of `hyper` are created: https://github.com/sympy/sympy/blob/a8a3a3b026cc55aa14010fc7cd7909806b6e116c/sympy/functions/special/hyper.py#L185-L187

`not(x0==0)` is usually written `x0 != 0`.
`return super(gamma, self)._eval_nseries(x, n, logx)`  should not contain gamma if this is for hyper. In fact, I think that `super()` should suffice as Python 2 need not be supported.

Otherwise, the series expansion looks correct (though I did not check carefully).
@jksuom I made the following changes but there are some tests failing locally.
```
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -220,6 +220,32 @@ def _eval_rewrite_as_Sum(self, ap, bq, z, **kwargs):
         return Piecewise((Sum(coeff * z**n / factorial(n), (n, 0, oo)),
                          self.convergence_statement), (self, True))

+    def _eval_nseries(self, x, n, logx):
+
+        from sympy.functions import factorial, RisingFactorial
+
+        x0 = self.args[2].limit(x, 0)
+        ap = self.args[0]
+        bq = self.args[1]
+
+        if x0 != 0:
+            return super()._eval_nseries(x, n, logx)
+
+        Rf1 = 1
+        Rf2 = 1
+        series = 0
+
+        for i in range(n):
+            for a in ap:
+                Rf1 *= RisingFactorial(a, i)
+
+            for b in bq:
+                Rf2 *= RisingFactorial(b, i)
+
+            series += ((Rf1 / Rf2) * (x ** i)) / factorial(i)
+
+        return series
+
     @property
     def argument(self):
         """ Argument of the hypergeometric function. """

```
```
________________________________________________________________________________
_______ sympy\functions\special\tests\test_elliptic_integrals.py:test_K ________
Traceback (most recent call last):
  File "c:\users\mendiratta\sympy\sympy\functions\special\tests\test_elliptic_integrals.py", line 41, in test_K
    25*pi*z**3/512 + 1225*pi*z**4/32768 + 3969*pi*z**5/131072 + O(z**6)
AssertionError
________________________________________________________________________________
_______ sympy\functions\special\tests\test_elliptic_integrals.py:test_E ________
Traceback (most recent call last):
  File "c:\users\mendiratta\sympy\sympy\functions\special\tests\test_elliptic_integrals.py", line 111, in test_E
    5*pi*z**3/512 - 175*pi*z**4/32768 - 441*pi*z**5/131072 + O(z**6)
AssertionError
________________________________________________________________________________
___________ sympy\functions\special\tests\test_hyper.py:test_limits ____________
Traceback (most recent call last):
  File "c:\users\mendiratta\sympy\sympy\functions\special\tests\test_hyper.py", line 347, in test_limits
    O(k**6) # issue 6350
AssertionError

 tests finished: 441 passed, 3 failed, 12 skipped, 10 expected to fail,
in 146.80 seconds
DO *NOT* COMMIT!
```
I can't understand why

Commit message: Implemented _eval_nseries() for Hyper

In reference to Issue #18193 the function was implemented. Example of
the working code is as follows:

In [1]: from sympy import hyper

In [2]: from sympy.abc import x, y

In [3]: hyper((1, 2),(1, 2), x)._eval_nseries(x, 7, y)
Out[3]: 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + x**6/720 + O(x**7)

In [4]: hyper((1, 2),(3, 4, 5), x)._eval_nseries(x, 7, y)
Out[4]: 1 + x/30 + x**2/1200 + x**3/63000 + x**4/4233600 + x**5/355622400 + x**6/36578304000 + O(x**7)

Patch to review:
diff --git a/sympy/functions/special/hyper.py b/sympy/functions/special/hyper.py
index ffda6bac2b0a..84d6153cfed6 100644
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -220,6 +220,33 @@ def _eval_rewrite_as_Sum(self, ap, bq, z, **kwargs):
         return Piecewise((Sum(coeff * z**n / factorial(n), (n, 0, oo)),
                          self.convergence_statement), (self, True))
 
+    def _eval_nseries(self, x, n, logx):
+
+        from sympy.functions import factorial, RisingFactorial
+        from sympy import Order, Add
+
+        x0 = self.args[2].limit(x, 0)
+        ap = self.args[0]
+        bq = self.args[1]
+
+        if x0 != 0:
+            return super()._eval_nseries(x, n, logx)
+
+        terms = []
+
+        for i in range(n):
+            Rf1 = 1
+            Rf2 = 1
+            for a in ap:
+                Rf1 *= RisingFactorial(a, i)
+
+            for b in bq:
+                Rf2 *= RisingFactorial(b, i)
+
+            terms.append(((Rf1 / Rf2) * (x ** i)) / factorial(i))
+
+        return (Add(*terms) + Order(x**n,x))
+
     @property
     def argument(self):
         """ Argument of the hypergeometric function. """


Code review comments:
diff_hunk: @@ -220,6 +220,33 @@ def _eval_rewrite_as_Sum(self, ap, bq, z, **kwargs):
         return Piecewise((Sum(coeff * z**n / factorial(n), (n, 0, oo)),
                          self.convergence_statement), (self, True))
 
+    def _eval_nseries(self, x, n, logx):
+
+        from sympy.functions import factorial, RisingFactorial
+        from sympy import Order, Add
+
+        x0 = self.args[2].limit(x, 0)
+        ap = self.args[0]
+        bq = self.args[1]
+
+        if x0 != 0:
+            return super()._eval_nseries(x, n, logx)
+
+        terms = []
+
+        for i in range(n):
+            Rf1 = 1
+            Rf2 = 1
+            for a in ap:
+                Rf1 *= RisingFactorial(a, i)
+
+            for b in bq:
+                Rf2 *= RisingFactorial(b, i)
+
+            terms.append(((Rf1 / Rf2) * (x ** i)) / factorial(i))
path: sympy/functions/special/hyper.py
text: @user1:
```suggestion
            terms.append(((Rf1/Rf2) * (arg**i)) / factorial(i))
```

diff_hunk: @@ -220,6 +220,33 @@ def _eval_rewrite_as_Sum(self, ap, bq, z, **kwargs):
         return Piecewise((Sum(coeff * z**n / factorial(n), (n, 0, oo)),
                          self.convergence_statement), (self, True))
 
+    def _eval_nseries(self, x, n, logx):
+
+        from sympy.functions import factorial, RisingFactorial
+        from sympy import Order, Add
+
+        x0 = self.args[2].limit(x, 0)
+        ap = self.args[0]
+        bq = self.args[1]
+
+        if x0 != 0:
+            return super()._eval_nseries(x, n, logx)
path: sympy/functions/special/hyper.py
text: @user1:
```suggestion
            return super()._eval_nseries(self, x, n, logx)
```

@user1:
I'm not sure of this, though. The docs say that if the second argument of `super` is omitted, then the method will be unbound but I don't know what will happen if both arguments are omitted. However, there is an example that seems to indicate that the returned method will be bound to `self`. Perhaps that should be tested.

@user1:
Tests seem to show that the method will be bound. So better leave this unchanged.

@author:
I think what you are saying is the way it was done in python 2 and in python 3 both the arguments can be left omitted.

diff_hunk: @@ -220,6 +220,33 @@ def _eval_rewrite_as_Sum(self, ap, bq, z, **kwargs):
         return Piecewise((Sum(coeff * z**n / factorial(n), (n, 0, oo)),
                          self.convergence_statement), (self, True))
 
+    def _eval_nseries(self, x, n, logx):
+
+        from sympy.functions import factorial, RisingFactorial
+        from sympy import Order, Add
+
+        x0 = self.args[2].limit(x, 0)
path: sympy/functions/special/hyper.py
text: @user1:
```suggestion
        arg = self.args[2]
        x0 = arg.limit(x, 0)
```