Repository: sympy/sympy

Issue:
limits of x*besselk(1,x) and x**3*besselk(1,x) as x -> 0 return implementation errors
hi, these should converge to one and zero respectively. was bummed to get an error while working out a derivation. perhaps there is a more general solution for taking the limit x**m*besselk(n,x) as x->0, but these two cases were important for what i was working on.

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* series
  * Fixed series expansions for bessely functions where terms from series were missing. ([#23844](https://github.com/sympy/sympy/pull/23844) by [@anutosh491](https://github.com/anutosh491))

  * Implemented _eval_as_leading_term method for besseli function. ([#23844](https://github.com/sympy/sympy/pull/23844) by [@anutosh491](https://github.com/anutosh491))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.12.

<details><summary>Click here to see the pull request description that was parsed.</summary>



    <!-- 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. -->


    #### Brief description of what is fixed or changed
    As of now , the pr implements leading term method for the besseli function and fixes series expansions for bessely function. The bessely function on master if neglecting few terms . For eg
    ```
    >>> # On master
    >>> bessely(1, x).series(x, 0, 3)  # first term -2/(pi*x) is missing
    x*(log(x)/pi - log(2)/pi - (1 - 2*EulerGamma)/(2*pi)) + O(x**3*log(x))
    >>>
    >>> bessely(2, x).series(x, 0, 3) # second term -1/pi  is missing
    -4/(pi*x**2) + x**2*(log(x)/(4*pi) - log(2)/(4*pi) - (3/2 - 2*EulerGamma)/(8*pi)) + O(x**3*log(x))
    >>>
    >>> bessely(3, x).series(x, 0, 3)   # Second term is wrong and also missing 3rd term x/(4*pi)
    -16/(pi*x**3) - 4/(pi*x) + O(x**3*(-log(x) + log(2)))
    >>>
    >>> # On branch
    >>> bessely(1, x).series(x, 0, 3)
    -2/(pi*x) + x*(log(x)/pi - log(2)/pi - (1 - 2*EulerGamma)/(2*pi)) + O(x**3*log(x))
    >>> bessely(2, x).series(x, 0, 3)
    -4/(pi*x**2) - 1/pi + x**2*(log(x)/(4*pi) - log(2)/(4*pi) - (3/2 - 2*EulerGamma)/(8*pi)) + O(x**3*log(x))
    >>> bessely(3, x).series(x, 0, 3)
    -16/(pi*x**3) - 2/(pi*x) - x/(4*pi) + O(x**3*(-log(x) + log(2)))
    ```

    #### 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.

    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
      * Fixed series expansions for bessely functions where terms from series were missing.
      * Implemented _eval_as_leading_term method for besseli function. 
    <!-- END RELEASE NOTES -->


</details><p>

github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
anutosh491: It seems this would also solve the series errors from #7224 ( the issue also needs fixes in the solvers module) but would solve the series issues reported there .
```
>>> # On master
>>> y1 = sqrt(x)*besseli(1,sqrt(8*x))
>>>
>>> limit(x*diff(y1,x,x)/y1,x,0) 
Traceback (most recent call last):
............. 
  File "C:\Users\91989\sympy\sympy\sympy\simplify\radsimp.py", line 1114, in fraction
    if term.is_commutative and (term.is_Pow or isinstance(term, exp)):
  File "C:\Users\91989\sympy\sympy\sympy\functions\elementary\exponential.py", line 214, in __instancecheck__
    return isinstance(instance, Pow) and instance.base is S.Exp1
RecursionError: maximum recursion depth exceeded while calling a Python object
>>>
>>>  # On branch
>>> limit(x*diff(y1,x,x)/y1,x,0)
2
```


anutosh491: After making the changes , there were 2 errors from `test_wester.py` which weren't working as expected. The latest commit fixes 1 of them . The error was
```
>>limit(log(x)/(log(x) + sin(x)), x, oo)
1
```
The expected ans is through heuristics . Originally as `calculate_series` was being used , calculate series correctly raises a `ValueError` for this from where the execution moves to the heuristic function, but now as are using leadterm , there were few places which needed improvements to fix this 

1) sign function in gruntz.py
```
def sign(e, x):
    if e.is_positive:
        return 1
    elif e.is_negative:
        return -1
    elif e.is_zero:
        return 0
    elif not e.has(x):
        from sympy.simplify import logcombine
        e = logcombine(e)
        return _sign(e)
```
The `_sign` function here is the sign function from `elementary/complexes.py` . So the limit involves an intermediate call to `sign(AccumBounds(-oo, oo), p)` which return `sign(AccumBounds(-oo, oo))` . This is the reason why `calculate_series` generates the `ValueError` but we currently haven't implemented the same in the sign class. I've addressed this through something like
```
+++ b/sympy/functions/elementary/complexes.py
@@ -328,6 +328,7 @@ def doit(self, **hints):

     @classmethod
     def eval(cls, arg):
+        from sympy.calculus.util import AccumBounds
@@ -375,6 +376,13 @@ def eval(cls, arg):
+        if isinstance(arg, AccumBounds):
+            if arg.max < S.Zero:
+                return S.NegativeOne
+            elif arg.min > S.Zero:
+                return S.One
+            else:
+                raise ValueError('sign(%s) cannot be determined' % (arg))
```
2) There is an intermediate call to the following  where (AccumBounds(-1, 1) is generated through sin(1/x).as_leading_term(x))
```
>>> 1/(-log(_w) + AccumBounds(-1, 1)).leadterm(w, logx = -p)
AccumBounds(-oo, oo)
```
The answer here isn't correct . Anything like `1/(p + AccumBounds(-1, 1))` might look better here This is eventually because of this 
```
>>> # Take a positive variable as _w defined in gruntz is a positive dummy
>>> y = Symbol('y', positive=True)
>>> y + AccumBounds(-1, 1)
AccumBounds(y - 1, y + 1)
>>> (_)**(-1)
AccumBounds(-oo, oo)
```
I guess this isn't correct cause y can be any positive number . I've addressed this through
```
+++ b/sympy/calculus/accumulationbounds.py
@@ -532,6 +535,8 @@ def __pow__(self, other):
                             return AccumBounds(self.max**other, oo)
                         if self.max.is_zero:
                             return AccumBounds(-oo, self.min**other)
+                        if self.has(Symbol):
+                            return AccumBounds(self.max**other, self.min**other)
```

3) Lastly there is an intermediate call to `gruntz(log(1/x) + sin(1/x), x, 0)` . As I've written above , as we're shifitng from `calculate_series` to `leadterm` ..... originally `calculate_series` would give `PoleError`  because of `sin(1/_w)._eval_nseries(_w)` but through leadterm we would eventually encounter an `AccumBounds(-1, 1)` term. If we don't take care of these cases , we'll get the following
```
>>> gruntz(log(1/x) + sin(1/x), x, 0)
Traceback (most recent call last):
NotImplementedError: Don't know how to calculate the mrv of 'AccumBounds(_p - 1, _p + 1)'
```
But like at the end it's not different then how mrv/mrv_leadterm should be operated for `p-1` or `p + 1` (obviously `gruntz(log(1/x) + 1, x, 0)` and `gruntz(log(1/x) - 1, x, 0)` do work) . Hence i've just added a block to fix this as of now 
```
+++ b/sympy/series/gruntz.py

@@ -442,7 +443,12 @@ def limitinf(e, x):
-    c0, e0 = mrv_leadterm(e, x)
+    if isinstance(e, AccumBounds):
+        if mrv_leadterm(e.min, x) != mrv_leadterm(e.max, x):
+            raise NotImplementedError
+        c0, e0 = mrv_leadterm(e.min, x)
+    else:
+        c0, e0 = mrv_leadterm(e, x)
```
This took some careful debugging to find out where exactly the errors turn up and I think I've managed to figure where they really turn up . The fixes I've added can surely be discussed . I expect a single error from `test_wester.py` to  show up now (simplification error)





jksuom: How does the addition of AccumBounds to `_sign` help with `gruntz.sign`?
anutosh491: > How does the addition of AccumBounds to `_sign` help with `gruntz.sign`?

Not a hard requirement, I put that just in case , thinking it might have a role to play for `limit(log(x)/(log(x) + sin(x)), x, oo)` but wasn't necessary . I've removed it in the latest commit.
jksuom: That AccumBounds code in `sign` looks correct but I suspect that it would not be of much practical significance.
jksuom: 2 and 3 look plausible to me. 
anutosh491: I guess we need to have a block to handle infinites better . As I was trying out asymptotic limits along directions other than just `+1/-1`..... I ended up trying these cases .
```
>>> # All outputs are wrong
>>>
>>> limit(besselj(0, x), x, -I*oo)
0
>>> limit(besselj(1, x), x, -I*oo)
0
>>> limit(besselj(2, x), x, -I*oo)
0
>>> limit(besselj(3, x), x, -I*oo)
0
```
As I am currently using the expression concerned with the lower plane, I think this code should be capable to handle these cases .
```
        if e.is_positive:
            return arg**nu/(2**nu*gamma(nu + 1))
        elif e.is_negative:
            return sqrt(2)*(-(-1)**(S(3)/4)*I**nu*exp(-I*pi*nu + I*z) + \
                    (-1)**(S(1)/4)*I**nu*exp(-I*z))/(2*sqrt(pi)*sqrt(z))
```
As of now we only make these adjustments for infinites other than `oo/-oo`
```
class Limit(Expr):

    def __new__(cls, e, z, z0, dir="+"):

        if z0 in (S.Infinity, S.ImaginaryUnit*S.Infinity):
            dir = "-"
        elif z0 in (S.NegativeInfinity, S.ImaginaryUnit*S.NegativeInfinity):
            dir = "+"
```
It's not that this is wrong , we should keep it but we need some more modifications in `Limits.doit()` block. We can solve this through
```
+++ b/sympy/series/limits.py
@@ -229,6 +229,13 @@ def doit(self, **hints):
         if z0 is S.ComplexInfinity:
             raise NotImplementedError("Limits at complex "
                                     "infinity are not implemented")
+
+        if z0.is_infinite:
+            cdir = sign(z0).simplify()
+            cdir = cdir/abs(cdir)
+            e = e.subs(z, cdir*z)
+            dir = "-"
+            z0 = S.Infinity
````
To get all correct answers as we want them
```
>>> limit(besselj(0, x), x, -I*oo)
oo
>>> limit(besselj(1, x), x, -I*oo)
-oo*I
>>> limit(besselj(2, x), x, -I*oo)
-oo
>>> limit(besselj(3, x), x, -I*oo)
oo*I
```
I'll check for failing tests if any and commit this patch

anutosh491: Limits for Shi function at `+oo/-oo` were giving wrong outputs due to a wrongly placed check in `Shi._eval_as_leading_term()`, the latest commit also fixes that .
```
>>> # On master
>>> limit(Shi(x), x, oo)   # expected oo
-I*pi/2
>>> limit(Shi(x), x, -oo)   # expected -oo
-I*pi/2
```
Reference from wolfram
<img width="78" alt="image" src="https://user-images.githubusercontent.com/87052487/184519026-fd2b900b-4717-4e73-8fcd-f8c2306d49fc.png">

anutosh491: @jksuom the latest commit is responsible for deprecating `calculate_series` function from `gruntz.py` . I haven't really dealt with deprecations before so would be glad if you could have a look at the changes made through the last commit . I've also replaced `calculate_series` with `series` at the couple of places where it was being used throughout the code base.
anutosh491: The failing test in **test / optional-dependencies (3.11.0-alpha - 3.11, true) (pull_request)** don't seem to be because of the code changes I've made today .
anutosh491: I also found that #6350 was pending because sympy lacked `_eval_nseries` implementation for `besseli` function . So currently we have the following on master
```
>>> k = Symbol('k')
>>> y = Symbol('y')
>>> integrate(exp(k*(y**3 - 3*y)), (y, 0, oo), conds='none')
-sqrt(3)*(k**(4/3)*exp(-2*I*pi/3)*gamma(-2/3)*gamma(-1/3)*hyper((1,), (4/3, 5/3), k**2) + k**(1/3)*besseli(-1/3, 2*k)*gamma(1/3)**2*gamma(2/3)**2 + k**(1/3)*exp(-I*pi/3)*besseli(1/3, 2*k)*gamma(-1/3)*gamma(1/3)*gamma(2/3)*gamma(4/3))*exp(2*I*pi/3)/(6*pi*k**(1/3))
>>>
>>> integrate(exp(k*(y**3 - 3*y)), (y, 0, oo), conds='none').series(k)
Traceback (most recent call last):
..........................................
  File "C:\Users\91989\sympy\sympy\sympy\core\function.py", line 736, in _eval_nseries
    raise PoleError("Cannot expand %s around 0" % (self))
sympy.core.function.PoleError: Cannot expand besseli(-1/3, 2*_x) around 0
```
But now that we have the implementation done , on the current branch , we get the expected series
```
>>> integrate(exp(k*(y**3 - 3*y)), (y, 0, oo), conds='none').series(k)
-(-1)**(2/3)*sqrt(3)*gamma(1/3)**2*gamma(2/3)/(6*pi*k**(1/3)) - sqrt(3)*k*gamma(-2/3)*gamma(-1/3)/(6*pi) - 3*sqrt(3)*k**3*gamma(-2/3)*gamma(-1/3)/(40*pi) - 27*sqrt(3)*k**5*gamma(-2/3)*gamma(-1/3)/(2240*pi) - (-1)**(1/3)*sqrt(3)*k**(1/3)*gamma(-1/3)*gamma(1/3)*gamma(2/3)/(6*pi) - (-1)**(2/3)*sqrt(3)*k**(5/3)*gamma(1/3)**2*gamma(2/3)/(4*pi) - (-1)**(1/3)*sqrt(3)*k**(7/3)*gamma(-1/3)*gamma(1/3)*gamma(2/3)/(8*pi) - 3*(-1)**(2/3)*sqrt(3)*k**(11/3)*gamma(1/3)**2*gamma(2/3)/(40*pi) - 3*(-1)**(1/3)*sqrt(3)*k**(13/3)*gamma(-1/3)*gamma(1/3)*gamma(2/3)/(112*pi) - 3*(-1)**(2/3)*sqrt(3)*k**(17/3)*gamma(1/3)**2*gamma(2/3)/(320*pi) + O(k**6)
>>>
>>> # The series expansion approximates the function pretty well
>>> integrate(exp(k*(y**3 - 3*y)), (y, 0, oo), conds='none').subs(k, 0.1).n()
1.14235423238222 - 1.14277316068123*I
>>> integrate(exp(k*(y**3 - 3*y)), (y, 0, oo), conds='none').series(k).removeO().subs(k, 0.1).n()
1.14235422817977 - 1.14277316937034*I
```

I'll add a test for the same to close the issue .
anutosh491: > This is a small change, though I'm not fully sure if there is something more, I need to take care of while handling these cases i.e., n = -k + 3, because I don't think we currently have any such test implemented in the code base .But looks straightforward so maybe I'll try committing after running some tests.

I've made the change . Two blocks mainly where the change was introduced , which are
```
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1677,16 +1677,19 @@ def _eval_nseries(self, x, n, logx, cdir=0):
             e = logcombine(e).cancel()

         if not (m.is_zero or e.is_number and e.is_real):
-            res = exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
-            if res is exp(e*log(b)):
-                return self
-            return res
+            if self == self._eval_as_leading_term(x, logx=logx, cdir=cdir):
+                res = exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
+                if res is exp(e*log(b)):
+                    return self
+                return res

@@ -1785,7 +1788,10 @@ def mul(d1, d2):

         if not (e.is_integer and e.is_positive and (e*d - n).is_nonpositive and
                 res == _mexpand(self)):
-            res += Order(x**n, x)
+            try:
+                res += Order(x**n, x)
+            except NotImplementedError:
+                return exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
         return res
```
I've explained why we would need to make the first change . The second change is simply because of what happens in this expression.
```
>>> x**y +O(x)
Traceback (most recent call last):
.............
NotImplementedError: Not sure of sign of 1 - y
```
In such cases we could return the series of the exponential form of self. (**NOTE** : Truly speaking even the exponential form shouldn't work because we are kind of doing the exact same operation , but there is a minute difference which leads to change in execution path of how both of these are dealt. The difference revolves around the following
```
>>> (x**y).leadterm(x)
(1, y)
>>> exp(y*log(x)).leadterm(x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3532, in leadterm
    l = self.as_leading_term(x, logx=logx, cdir=cdir)
  File "C:\Users\91989\sympy\sympy\sympy\core\cache.py", line 70, in wrapper
    retval = cfunc(*args, **kwargs)
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3495, in as_leading_term
    obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
  File "C:\Users\91989\sympy\sympy\sympy\functions\elementary\exponential.py", line 551, in _eval_as_leading_term
    raise PoleError("Cannot expand %s around 0" % (self))
sympy.core.function.PoleError: Cannot expand exp(y*log(x)) around 0
```
)

Hence something like `(sin(x)**y).nseries(x)` can be answered , just as we have been doing on master
```
>>> ((sin(x))**y).nseries(x, n=1, logx=Symbol('logx'))
exp(logx*y) + O(x)
>>> ((sin(x))**y).nseries(x, n=6, logx=Symbol('logx'))
exp(logx*y) - x**2*y*exp(logx*y)/6 - x**4*y*exp(logx*y)/180 + x**4*y**2*exp(logx*y)/72 + O(x**6)




anutosh491: Hello @Parcly-Taxel , would you help me fix the failing test here . It came into picture because #23985 was merged some days back. The Pr involved changes in `_eval_rewrite_as_tractable` function of harmonic class.
```
    def _eval_rewrite_as_tractable(self, n, m=1, limitvar=None, **kwargs):
        from sympy.functions.special.zeta_functions import zeta
        from sympy.functions.special.gamma_functions import polygamma
        pg = self.rewrite(polygamma)
        if not isinstance(pg, harmonic):
            return pg.rewrite("tractable", deep=True)
        return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
```
Now for limits like these....
```
>>> x = Symbol('x')
>>> y = Symbol('y', positive=True)
>>> limit(harmonic(x, y), x, oo)
harmonic(oo, y)
>>> limit(harmonic(x, y + 1), x, oo)
zeta(y + 1)
```
... On master the results are being returned through `heuristics` ve because `calculate_series` method from `gruntz` raises `PoleError` for these
```
>>> gruntz(harmonic(x, y), x, oo)
Traceback (most recent call last):
.........................
sympy.core.function.PoleError:
Asymptotic expansion of zeta around [y, oo] is not implemented.
>>>
>>> gruntz(harmonic(x, y + 1), x, oo)
Traceback (most recent call last):
..............
sympy.core.function.PoleError:
Asymptotic expansion of zeta around [y + 1, oo] is not implemented.
```

But this pr is working towards deprecation of the calculate_series method as it is buggy quite sometimes and inefficient too and proposed moving towards `leadterm/as_leading_term` methods which should be good enough to do the job.

Hence now that your Pull request changed `_eval_rewrite_as_tracatable` method for the harmonic class , there is a failing tests on the pr . It goes as follows
```
>>> y = Symbol('y', positive=True)
>>> limit(harmonic(x, y), x, oo)
zeta(y)
```

This is because now that we don't use `calculate_series`, we won't be getting the `PoleError` and we would be getting the answer out of the `gruntz` call rather than `heursitics` . Now in gruntz.py , we have encounter the following code execution
```
# initially e  is  harmonic(_p, y)
# line 446 is executed

e = e.rewrite('tractable', deep=True, limitvar=x)

# e is changed to zeta(y) - zeta(y, _p + 1)
# line 454 is executed
# this calculates mrv_leadterm of zeta(y) - zeta(y, _p + 1)
# the answer for this would be (zeta(y), 0) because of the following expression (which is correct) 

>>> zeta(y) - zeta(y, 1 + 1/x)
zeta(y) - zeta(y, 1 + 1/x)
>>> _.as_leading_term(x)
zeta(y)

# Now line 469 is executed and zeta(y) is returned as answer
```
I have been looking in this since yesterday and I can see a couple of ways which could fix this , but I am not sure if it would be best done that way . Hence any input from your side would be highly appreciated !
 


anutosh491: Using the rewrite function as tractable functions which we have currently might not always work .
```
    def _eval_rewrite_as_tractable(self, n, m=1, limitvar=None, **kwargs):
        if not isinstance(pg, harmonic):
            return pg.rewrite("tractable", deep=True)
        return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
```
As then we would always have
```
>>> x = Symbol('x')
>>> y = Symbol('y', positive=True)
>>> z = Symbol('z', negative=True)
>>> limit(harmonic(x, y + 1), x, oo)
zeta(y + 1)
>>> limit(harmonic(x, y), x, oo)
zeta(y)
>>> limit(harmonic(x, z - 1), x, oo)
zeta(z - 1)
>>> limit(harmonic(x,  z), x, oo)
zeta(z)
```
Becuase this result would always depend on 
`(zeta(f(y)) - zeta(f(y)), 1 + 1/x)).as_leading_term(x)` which is always `zeta(f(y))` . So could you maybe have a look into the `_eval_rewrite_as_tractable` function which I am pointing out @Parcly-Taxel 
anutosh491: Okay I realize there is a small problem with how the code flow is executed . What's wrong with the general code is the following
```
>>> x, y = symbols('x y')
>>> zeta(y, 1 + 1/x).as_leading_term(x)
zeta(y, 1/x)
```
Now this is only completely true when `y > 1` , below that we can't really comment unless we know more about values y can take.
For eg
```
>>> zeta(3, 1 + 1/x).as_leading_term(x)
zeta(3, 1/x)
>>> zeta(4, 1 + 1/x).as_leading_term(x)
zeta(4, 1/x)
>>> zeta(1, 1 + 1/x).as_leading_term(x)
zoo
>>> zeta(0, 1 + 1/x).as_leading_term(x)
-1/x
>>> zeta(-1, 1 + 1/x).as_leading_term(x)
-1/(2*x**2)
```
Hence we can't really allow the general case i.e. `zeta(y, 1 + 1/x).as_leading_term(x)`  to return `zeta(y, 1/x)` everytime unless and untill there are any assumptions on `y` (either positive , negative or something) .
How I plan to fix this 

1) First fix `_eval_rewrite_as_tractable`
As we know the relation between harmonic and the zeta function is 
$$H_{n,m}=\begin{cases}\psi(n+1)+\gamma&m=1\\\\ \zeta(m)-\zeta(m,n+1)&m\ne1\end{cases}$$

We need to use something like the following
```
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -988,7 +988,11 @@ def _eval_rewrite_as_tractable(self, n, m=1, limitvar=None, **kwargs):
         pg = self.rewrite(polygamma)
         if not isinstance(pg, harmonic):
             return pg.rewrite("tractable", deep=True)
-        return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
+        arg = m - S.One
+        if arg.is_positive or arg.is_negative:
+            return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
```

2) We need to take care of that general case . For now I think we can introduce a maybe not so full fledged _eval_as_leading_term method for zeta function . Something like the following
```
--- a/sympy/functions/special/zeta_functions.py
+++ b/sympy/functions/special/zeta_functions.py
@@ -561,6 +561,20 @@ def fdiff(self, argindex=1):
             return -s*zeta(s + 1, a)
         else:
             raise ArgumentIndexError
+
+    def _eval_as_leading_term(self, x, logx=None, cdir=0):
+
+        s, a = self.args
+        try:
+            c, e = a.leadterm(x)
+        except NotImplementedError:
+            return self
+
+        if e.is_negative and not s.is_positive:
+            raise NotImplementedError
+
+        return super(zeta, self)._eval_as_leading_term(x, logx, cdir)
```
This is will give us a range of expected asnwers
```
>>> x = Symbol('x')
>>> y = Symbol('y', positive=True)
>>> z = Symbol('z', negative=True)
>>> limit(harmonic(x, y), x, oo)
harmonic(oo, y)
>>> limit(harmonic(x, y + 1), x, oo)
zeta(y + 1)
>>> limit(harmonic(x, y + 2), x, oo)
zeta(y + 2)
>>> limit(harmonic(x, y + 2), x, oo)
zeta(y + 2)
>>> limit(harmonic(x, y - 1), x, oo)
harmonic(oo, y - 1)
>>> limit(harmonic(x, z), x, oo)
Limit(harmonic(x, z), x, oo, dir='-')
>>> limit(harmonic(x, z + 1), x, oo)
oo
>>> limit(harmonic(x, z + 2), x, oo)
harmonic(oo, z + 2)
>>> limit(harmonic(x, z + 3), x, oo)
harmonic(oo, z + 3)
>>> limit(harmonic(x, z - 1), x, oo)
Limit(harmonic(x, z - 1), x, oo, dir='-')
>>> limit(harmonic(x, z - 2), x, oo)
Limit(harmonic(x, z - 2), x, oo, dir='-')
```



Patch:
diff --git a/sympy/calculus/accumulationbounds.py b/sympy/calculus/accumulationbounds.py
index 93fa7386e7..44f2272a0d 100644
--- a/sympy/calculus/accumulationbounds.py
+++ b/sympy/calculus/accumulationbounds.py
@@ -523,7 +523,7 @@ def __pow__(self, other):
                             return AccumBounds(self.max**other, oo)
                         if self.max.is_zero:
                             return AccumBounds(self.min**other, oo)
-                        return AccumBounds(0, oo)
+                        return (1/self)**(-other)
                     return AccumBounds(
                         S.Zero, Max(self.min**other, self.max**other))
                 elif other % 2 == 1:
@@ -532,7 +532,7 @@ def __pow__(self, other):
                             return AccumBounds(self.max**other, oo)
                         if self.max.is_zero:
                             return AccumBounds(-oo, self.min**other)
-                        return AccumBounds(-oo, oo)
+                        return (1/self)**(-other)
                     return AccumBounds(self.min**other, self.max**other)
 
             # non-integer exponent
diff --git a/sympy/core/add.py b/sympy/core/add.py
index ba6bf85fd9..aebc468ecd 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -1009,7 +1009,7 @@ def as_real_imag(self, deep=True, **hints):
         return (self.func(*re_part), self.func(*im_part))
 
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
-        from sympy.core.symbol import Dummy
+        from sympy.core.symbol import Dummy, Symbol
         from sympy.series.order import Order
         from sympy.functions.elementary.exponential import log
         from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
@@ -1067,7 +1067,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             try:
                 n0 = min.getn()
             except NotImplementedError:
-                n0 = 1
+                n0 = S.One
+            if n0.has(Symbol):
+                n0 = S.One
             res = Order(1)
             incr = S.One
             while res.is_Order:
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index ba9ad0966d..40a221bfd5 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -3425,6 +3425,14 @@ def compute_leading_term(self, x, logx=None):
         as_leading_term is only allowed for results of .series()
         This is a wrapper to compute a series first.
         """
+        from sympy.utilities.exceptions import SymPyDeprecationWarning
+
+        SymPyDeprecationWarning(
+            feature="compute_leading_term",
+            useinstead="as_leading_term",
+            issue=21843,
+            deprecated_since_version="1.12"
+        ).warn()
 
         from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
         if self.has(Piecewise):
@@ -3434,17 +3442,22 @@ def compute_leading_term(self, x, logx=None):
         if self.removeO() == 0:
             return self
 
-        from sympy.series.gruntz import calculate_series
+        from .symbol import Dummy
+        from sympy.functions.elementary.exponential import log
+        from sympy.series.order import Order
+
+        _logx = logx
+        logx = Dummy('logx') if logx is None else logx
+        res = Order(1)
+        incr = S.One
+        while res.is_Order:
+            res = expr._eval_nseries(x, n=1+incr, logx=logx).cancel().powsimp().trigsimp()
+            incr *= 2
 
-        if logx is None:
-            from .symbol import Dummy
-            from sympy.functions.elementary.exponential import log
-            d = Dummy('logx')
-            s = calculate_series(expr, x, d).subs(d, log(x))
-        else:
-            s = calculate_series(expr, x, logx)
+        if _logx is None:
+            res = res.subs(logx, log(x))
 
-        return s.as_leading_term(x)
+        return res.as_leading_term(x)
 
     @cacheit
     def as_leading_term(self, *symbols, logx=None, cdir=0):
diff --git a/sympy/core/function.py b/sympy/core/function.py
index b68c865547..019b0d17e7 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -687,9 +687,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):
         args0 = [t.limit(x, 0) for t in args]
         if any(t.is_finite is False for t in args0):
             from .numbers import oo, zoo, nan
-            # XXX could use t.as_leading_term(x) here but it's a little
-            # slower
-            a = [t.compute_leading_term(x, logx=logx) for t in args]
+            a = [t.as_leading_term(x, logx=logx) for t in args]
             a0 = [t.limit(x, 0) for t in a]
             if any(t.has(oo, -oo, zoo, nan) for t in a0):
                 return self._eval_aseries(n, args0, x, logx)
diff --git a/sympy/core/power.py b/sympy/core/power.py
index b6f2ba4fe8..95f7c7ff7a 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1633,6 +1633,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):
         from sympy.functions.elementary.exponential import exp, log
         from sympy.series.limits import limit
         from sympy.series.order import Order
+        from sympy.core.sympify import sympify
         if self.base is S.Exp1:
             e_series = self.exp.nseries(x, n=n, logx=logx)
             if e_series.is_Order:
@@ -1686,16 +1687,19 @@ def _eval_nseries(self, x, n, logx, cdir=0):
             e = logcombine(e).cancel()
 
         if not (m.is_zero or e.is_number and e.is_real):
-            res = exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
-            if res is exp(e*log(b)):
-                return self
-            return res
+            if self == self._eval_as_leading_term(x, logx=logx, cdir=cdir):
+                res = exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
+                if res == exp(e*log(b)):
+                    return self
+                return res
 
         f = b.as_leading_term(x, logx=logx)
         g = (b/f - S.One).cancel(expand=False)
         if not m.is_number:
             raise NotImplementedError()
         maxpow = n - m*e
+        if maxpow.has(Symbol):
+            maxpow = sympify(n)
 
         if maxpow.is_negative:
             return Order(x**(m*e), x)
@@ -1794,7 +1798,10 @@ def mul(d1, d2):
 
         if not (e.is_integer and e.is_positive and (e*d - n).is_nonpositive and
                 res == _mexpand(self)):
-            res += Order(x**n, x)
+            try:
+                res += Order(x**n, x)
+            except NotImplementedError:
+                return exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
         return res
 
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
diff --git a/sympy/functions/combinatorial/numbers.py b/sympy/functions/combinatorial/numbers.py
index f8dea4642c..a58e0cae86 100644
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -984,7 +984,9 @@ def _eval_rewrite_as_tractable(self, n, m=1, limitvar=None, **kwargs):
         pg = self.rewrite(polygamma)
         if not isinstance(pg, harmonic):
             return pg.rewrite("tractable", deep=True)
-        return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
+        arg = m - S.One
+        if arg.is_nonzero:
+            return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
 
     def _eval_evalf(self, prec):
         if not all(x.is_number for x in self.args):
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 1302454473..65c406a632 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -1267,7 +1267,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):
             from sympy.functions.elementary.integers import ceiling
             from sympy.series.order import Order
             arg = self.args[0].nseries(x, n=n, logx=logx)
-            lt = arg.compute_leading_term(x, logx=logx)
+            lt = arg.as_leading_term(x, logx=logx)
             lte = 1
             if lt.is_Pow:
                 lte = lt.exp
diff --git a/sympy/functions/special/bessel.py b/sympy/functions/special/bessel.py
index 3c2a07a76c..ba0b57fcb1 100644
--- a/sympy/functions/special/bessel.py
+++ b/sympy/functions/special/bessel.py
@@ -223,11 +223,24 @@ def _eval_rewrite_as_jn(self, nu, z, **kwargs):
 
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
         nu, z = self.args
-        arg = z.as_leading_term(x)
-        if x in arg.free_symbols:
+        try:
+            arg = z.as_leading_term(x)
+        except NotImplementedError:
+            return self
+        c, e = arg.as_coeff_exponent(x)
+
+        if e.is_positive:
             return arg**nu/(2**nu*gamma(nu + 1))
-        else:
-            return self.func(nu, z.subs(x, 0))
+        elif e.is_negative:
+            cdir = 1 if cdir == 0 else cdir
+            sign = c*cdir**e
+            if not sign.is_negative:
+                # Refer Abramowitz and Stegun 1965, p. 364 for more information on
+                # asymptotic approximation of besselj function.
+                return sqrt(2)*cos(z - pi*(2*nu + 1)/4)/sqrt(pi*z)
+            return self
+
+        return super(besselj, self)._eval_as_leading_term(x, logx, cdir)
 
     def _eval_is_extended_real(self):
         nu, z = self.args
@@ -235,6 +248,8 @@ def _eval_is_extended_real(self):
             return True
 
     def _eval_nseries(self, x, n, logx, cdir=0):
+        # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
+        # for more information on nseries expansion of besselj function.
         from sympy.series.order import Order
         nu, z = self.args
 
@@ -319,6 +334,10 @@ def eval(cls, nu, z):
                 return S.NaN
         if z in (S.Infinity, S.NegativeInfinity):
             return S.Zero
+        if z == S.ImaginaryUnit*S.Infinity:
+            return exp(I*pi*(nu + 1)/2) * S.Infinity
+        if z == S.ImaginaryUnit*S.NegativeInfinity:
+            return exp(-I*pi*(nu + 1)/2) * S.Infinity
 
         if nu.is_integer:
             if nu.could_extract_minus_sign():
@@ -338,14 +357,28 @@ def _eval_rewrite_as_yn(self, nu, z, **kwargs):
 
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
         nu, z = self.args
-        term_one = ((2/pi)*log(z/2)*besselj(nu, z))
-        term_two = -(z/2)**(-nu)*factorial(nu - 1)/pi if (nu).is_positive else S.Zero
-        term_three = -(z/2)**nu/(pi*factorial(nu))*(digamma(nu + 1) - S.EulerGamma)
-        arg = Add(*[term_one, term_two, term_three]).as_leading_term(x, logx=logx)
-        if (x in arg.free_symbols or logx in arg.free_symbols):
+        try:
+            arg = z.as_leading_term(x)
+        except NotImplementedError:
+            return self
+        c, e = arg.as_coeff_exponent(x)
+
+        if e.is_positive:
+            term_one = ((2/pi)*log(z/2)*besselj(nu, z))
+            term_two = -(z/2)**(-nu)*factorial(nu - 1)/pi if (nu).is_positive else S.Zero
+            term_three = -(z/2)**nu/(pi*factorial(nu))*(digamma(nu + 1) - S.EulerGamma)
+            arg = Add(*[term_one, term_two, term_three]).as_leading_term(x, logx=logx)
             return arg
-        else:
-            return self.func(nu, z.subs(x, 0).cancel())
+        elif e.is_negative:
+            cdir = 1 if cdir == 0 else cdir
+            sign = c*cdir**e
+            if not sign.is_negative:
+                # Refer Abramowitz and Stegun 1965, p. 364 for more information on
+                # asymptotic approximation of bessely function.
+                return sqrt(2)*(-sin(pi*nu/2 - z + pi/4) + 3*cos(pi*nu/2 - z + pi/4)/(8*z))*sqrt(1/z)/sqrt(pi)
+            return self
+
+        return super(bessely, self)._eval_as_leading_term(x, logx, cdir)
 
     def _eval_is_extended_real(self):
         nu, z = self.args
@@ -353,6 +386,8 @@ def _eval_is_extended_real(self):
             return True
 
     def _eval_nseries(self, x, n, logx, cdir=0):
+        # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/0008/
+        # for more information on nseries expansion of bessely function.
         from sympy.series.order import Order
         nu, z = self.args
 
@@ -375,11 +410,15 @@ def _eval_nseries(self, x, n, logx, cdir=0):
                 return o
             t = (_mexpand(r**2) + o).removeO()
 
-            if nu > S.One:
+            if nu > S.Zero:
                 term = r**(-nu)*factorial(nu - 1)/pi
                 b.append(term)
-                for k in range(1, nu - 1):
-                    term *= t*(nu - k - 1)/k
+                for k in range(1, nu):
+                    denom = (nu - k)*k
+                    if denom == S.Zero:
+                        term *= t/k
+                    else:
+                        term *= t/denom
                     term = (_mexpand(term) + o).removeO()
                     b.append(term)
 
@@ -452,6 +491,10 @@ def eval(cls, nu, z):
                 return S.NaN
         if im(z) in (S.Infinity, S.NegativeInfinity):
             return S.Zero
+        if z is S.Infinity:
+            return S.Infinity
+        if z is S.NegativeInfinity:
+            return (-1)**nu*S.Infinity
 
         if z.could_extract_minus_sign():
             return (z)**nu*(-z)**(-nu)*besseli(nu, -z)
@@ -491,6 +534,58 @@ def _eval_is_extended_real(self):
         if nu.is_integer and z.is_extended_real:
             return True
 
+    def _eval_as_leading_term(self, x, logx=None, cdir=0):
+        nu, z = self.args
+        try:
+            arg = z.as_leading_term(x)
+        except NotImplementedError:
+            return self
+        c, e = arg.as_coeff_exponent(x)
+
+        if e.is_positive:
+            return arg**nu/(2**nu*gamma(nu + 1))
+        elif e.is_negative:
+            cdir = 1 if cdir == 0 else cdir
+            sign = c*cdir**e
+            if not sign.is_negative:
+                # Refer Abramowitz and Stegun 1965, p. 377 for more information on
+                # asymptotic approximation of besseli function.
+                return exp(z)/sqrt(2*pi*z)
+            return self
+
+        return super(besseli, self)._eval_as_leading_term(x, logx, cdir)
+
+    def _eval_nseries(self, x, n, logx, cdir=0):
+        # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselI/06/01/04/01/01/0003/
+        # for more information on nseries expansion of besseli function.
+        from sympy.series.order import Order
+        nu, z = self.args
+
+        # In case of powers less than 1, number of terms need to be computed
+        # separately to avoid repeated callings of _eval_nseries with wrong n
+        try:
+            _, exp = z.leadterm(x)
+        except (ValueError, NotImplementedError):
+            return self
+
+        if exp.is_positive:
+            newn = ceiling(n/exp)
+            o = Order(x**n, x)
+            r = (z/2)._eval_nseries(x, n, logx, cdir).removeO()
+            if r is S.Zero:
+                return o
+            t = (_mexpand(r**2) + o).removeO()
+
+            term = r**nu/gamma(nu + 1)
+            s = [term]
+            for k in range(1, (newn + 1)//2):
+                term *= t/(k*(nu + k))
+                term = (_mexpand(term) + o).removeO()
+                s.append(term)
+            return Add(*s) + o
+
+        return super(besseli, self)._eval_nseries(x, n, logx, cdir)
+
 
 class besselk(BesselBase):
     r"""
@@ -573,6 +668,76 @@ def _eval_is_extended_real(self):
         if nu.is_integer and z.is_positive:
             return True
 
+    def _eval_as_leading_term(self, x, logx=None, cdir=0):
+        nu, z = self.args
+        try:
+            arg = z.as_leading_term(x)
+        except NotImplementedError:
+            return self
+        _, e = arg.as_coeff_exponent(x)
+
+        if e.is_positive:
+            term_one = ((-1)**(nu -1)*log(z/2)*besseli(nu, z))
+            term_two = (z/2)**(-nu)*factorial(nu - 1)/2 if (nu).is_positive else S.Zero
+            term_three = (-1)**nu*(z/2)**nu/(2*factorial(nu))*(digamma(nu + 1) - S.EulerGamma)
+            arg = Add(*[term_one, term_two, term_three]).as_leading_term(x, logx=logx)
+            return arg
+        elif e.is_negative:
+            # Refer Abramowitz and Stegun 1965, p. 378 for more information on
+            # asymptotic approximation of besselk function.
+            return sqrt(pi)*exp(-z)/sqrt(2*z)
+
+        return super(besselk, self)._eval_as_leading_term(x, logx, cdir)
+
+    def _eval_nseries(self, x, n, logx, cdir=0):
+        # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselK/06/01/04/01/02/0008/
+        # for more information on nseries expansion of besselk function.
+        from sympy.series.order import Order
+        nu, z = self.args
+
+        # In case of powers less than 1, number of terms need to be computed
+        # separately to avoid repeated callings of _eval_nseries with wrong n
+        try:
+            _, exp = z.leadterm(x)
+        except (ValueError, NotImplementedError):
+            return self
+
+        if exp.is_positive and nu.is_integer:
+            newn = ceiling(n/exp)
+            bn = besseli(nu, z)
+            a = ((-1)**(nu - 1)*log(z/2)*bn)._eval_nseries(x, n, logx, cdir)
+
+            b, c = [], []
+            o = Order(x**n, x)
+            r = (z/2)._eval_nseries(x, n, logx, cdir).removeO()
+            if r is S.Zero:
+                return o
+            t = (_mexpand(r**2) + o).removeO()
+
+            if nu > S.Zero:
+                term = r**(-nu)*factorial(nu - 1)/2
+                b.append(term)
+                for k in range(1, nu):
+                    denom = (k - nu)*k
+                    if denom == S.Zero:
+                        term *= t/k
+                    else:
+                        term *= t/denom
+                    term = (_mexpand(term) + o).removeO()
+                    b.append(term)
+
+            p = r**nu*(-1)**nu/(2*factorial(nu))
+            term = p*(digamma(nu + 1) - S.EulerGamma)
+            c.append(term)
+            for k in range(1, (newn + 1)//2):
+                p *= t/(k*(k + nu))
+                p = (_mexpand(p) + o).removeO()
+                term = p*(digamma(k + nu + 1) + digamma(k + 1))
+                c.append(term)
+            return a + Add(*b) + Add(*c) # Order term comes from a
+
+        return super(besselk, self)._eval_nseries(x, n, logx, cdir)
+
 
 class hankel1(BesselBase):
     r"""
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index b3b4572ff8..4ab4247df4 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -2021,7 +2021,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if arg0 is S.NaN:
             arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
         if arg0.is_zero:
-            return S.EulerGamma
+            c, e = arg.as_coeff_exponent(x)
+            logx = log(x) if logx is None else logx
+            return log(c) + e*logx + S.EulerGamma
         elif arg0.is_finite:
             return self.func(arg0)
         else:
@@ -2146,8 +2148,6 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             return arg
         elif not arg0.is_infinite:
             return self.func(arg0)
-        elif arg0.is_infinite:
-            return -pi*S.ImaginaryUnit/2
         else:
             return self
 
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 67f5af8156..2b1dd1e7db 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -850,7 +850,8 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         n, z = [a.as_leading_term(x) for a in self.args]
         o = Order(z, x)
         if n == 0 and o.contains(1/x):
-            return o.getn() * log(x)
+            logx = log(x) if logx is None else logx
+            return o.getn() * logx
         else:
             return self.func(n, z)
 
diff --git a/sympy/functions/special/zeta_functions.py b/sympy/functions/special/zeta_functions.py
index 36016ca579..6c672e48fb 100644
--- a/sympy/functions/special/zeta_functions.py
+++ b/sympy/functions/special/zeta_functions.py
@@ -565,6 +565,22 @@ def fdiff(self, argindex=1):
         else:
             raise ArgumentIndexError
 
+    def _eval_as_leading_term(self, x, logx=None, cdir=0):
+        if len(self.args) == 2:
+            s, a = self.args
+        else:
+            s, a = self.args + (S.One,)
+
+        try:
+            c, e = a.leadterm(x)
+        except NotImplementedError:
+            return self
+
+        if e.is_negative and not s.is_positive:
+            raise NotImplementedError
+
+        return super(zeta, self)._eval_as_leading_term(x, logx, cdir)
+
 
 class dirichlet_eta(Function):
     r"""
diff --git a/sympy/series/gruntz.py b/sympy/series/gruntz.py
index 0e19e56b88..2065b9c528 100644
--- a/sympy/series/gruntz.py
+++ b/sympy/series/gruntz.py
@@ -118,7 +118,7 @@
 """
 from functools import reduce
 
-from sympy.core import Basic, S, Mul, PoleError
+from sympy.core import Basic, S, Mul, PoleError, expand_mul
 from sympy.core.cache import cacheit
 from sympy.core.numbers import ilcm, I, oo
 from sympy.core.symbol import Dummy, Wild
@@ -126,6 +126,7 @@
 
 from sympy.functions import log, exp, sign as _sign
 from sympy.series.order import Order
+from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.misc import debug_decorator as debug
 from sympy.utilities.timeutils import timethis
 
@@ -423,21 +424,15 @@ def sign(e, x):
 @debug
 @timeit
 @cacheit
-def limitinf(e, x, leadsimp=False):
-    """Limit e(x) for x-> oo.
-
-    Explanation
-    ===========
-
-    If ``leadsimp`` is True, an attempt is made to simplify the leading
-    term of the series expansion of ``e``. That may succeed even if
-    ``e`` cannot be simplified.
-    """
+def limitinf(e, x):
+    """Limit e(x) for x-> oo."""
     # rewrite e in terms of tractable functions only
 
+    old = e
     if not e.has(x):
         return e  # e is a constant
     from sympy.simplify.powsimp import powdenest
+    from sympy.calculus.util import AccumBounds
     if e.has(Order):
         e = e.expand().removeO()
     if not x.is_positive or x.is_integer:
@@ -449,7 +444,12 @@ def limitinf(e, x, leadsimp=False):
         x = p
     e = e.rewrite('tractable', deep=True, limitvar=x)
     e = powdenest(e)
-    c0, e0 = mrv_leadterm(e, x)
+    if isinstance(e, AccumBounds):
+        if mrv_leadterm(e.min, x) != mrv_leadterm(e.max, x):
+            raise NotImplementedError
+        c0, e0 = mrv_leadterm(e.min, x)
+    else:
+        c0, e0 = mrv_leadterm(e, x)
     sig = sign(e0, x)
     if sig == 1:
         return S.Zero  # e0>0: lim f = 0
@@ -462,9 +462,9 @@ def limitinf(e, x, leadsimp=False):
             raise ValueError("Leading term should not be 0")
         return s*oo
     elif sig == 0:
-        if leadsimp:
-            c0 = c0.simplify()
-        return limitinf(c0, x, leadsimp)  # e0=0: lim f = lim c0
+        if c0 == old:
+            c0 = c0.cancel()
+        return limitinf(c0, x)  # e0=0: lim f = lim c0
     else:
         raise ValueError("{} could not be evaluated".format(sig))
 
@@ -489,6 +489,14 @@ def calculate_series(e, x, logx=None):
 
     This is a place that fails most often, so it is in its own function.
     """
+
+    SymPyDeprecationWarning(
+        feature="calculate_series",
+        useinstead="series() with suitable n, or as_leading_term",
+        issue=21838,
+        deprecated_since_version="1.12"
+    ).warn()
+
     from sympy.simplify.powsimp import powdenest
 
     for t in e.lseries(x, logx=logx):
@@ -544,18 +552,24 @@ def mrv_leadterm(e, x):
     #
     w = Dummy("w", positive=True)
     f, logw = rewrite(exps, Omega, x, w)
-    series = calculate_series(f, w, logx=logw)
     try:
-        lt = series.leadterm(w, logx=logw)
-    except (ValueError, PoleError):
-        lt = f.as_coeff_exponent(w)
-        # as_coeff_exponent won't always split in required form. It may simply
-        # return (f, 0) when a better form may be obtained. Example (-x)**(-pi)
-        # can be written as (-1**(-pi), -pi) which as_coeff_exponent does not return
-        if lt[0].has(w):
-            base = f.as_base_exp()[0].as_coeff_exponent(w)
-            ex = f.as_base_exp()[1]
-            lt = (base[0]**ex, base[1]*ex)
+        lt = f.leadterm(w, logx=logw)
+    except (NotImplementedError, PoleError, ValueError):
+        n0 = 1
+        _series = Order(1)
+        incr = S.One
+        while _series.is_Order:
+            _series = f._eval_nseries(w, n=n0+incr, logx=logw)
+            incr *= 2
+        series = _series.expand().removeO()
+        try:
+            lt = series.leadterm(w, logx=logw)
+        except (NotImplementedError, PoleError, ValueError):
+            lt = f.as_coeff_exponent(w)
+            if lt[0].has(w):
+                base = f.as_base_exp()[0].as_coeff_exponent(w)
+                ex = f.as_base_exp()[1]
+                lt = (base[0]**ex, base[1]*ex)
     return (lt[0].subs(log(w), logw), lt[1])
 
 
@@ -609,6 +623,8 @@ def rewrite(e, Omega, x, wsym):
     Returns the rewritten e in terms of w and log(w). See test_rewrite1()
     for examples and correct results.
     """
+
+    from sympy import AccumBounds
     if not isinstance(Omega, SubsSet):
         raise TypeError("Omega should be an instance of SubsSet")
     if len(Omega) == 0:
@@ -627,7 +643,7 @@ def rewrite(e, Omega, x, wsym):
     # g is going to be the "w" - the simplest one in the mrv set
     for g, _ in Omega:
         sig = sign(g.exp, x)
-        if sig != 1 and sig != -1:
+        if sig != 1 and sig != -1 and not sig.has(AccumBounds):
             raise NotImplementedError('Result depends on the sign of %s' % sig)
     if sig == 1:
         wsym = 1/wsym  # if g goes to oo, substitute 1/w
@@ -670,6 +686,12 @@ def rewrite(e, Omega, x, wsym):
     f = f.subs({wsym: wsym**exponent})
     logw /= exponent
 
+    # bottom_up function is required for a specific case - when f is
+    # -exp(p/(p + 1)) + exp(-p**2/(p + 1) + p). No current simplification
+    # methods reduce this to 0 while not expanding polynomials.
+    f = bottom_up(f, lambda w: getattr(w, 'normal', lambda: w)())
+    f = expand_mul(f)
+
     return f, logw
 
 
@@ -707,10 +729,7 @@ def gruntz(e, z, z0, dir="+"):
         else:
             raise NotImplementedError("dir must be '+' or '-'")
 
-    try:
-        r = limitinf(e0, z)
-    except ValueError:
-        r = limitinf(e0, z, leadsimp=True)
+    r = limitinf(e0, z)
 
     # This is a bit of a heuristic for nice results... we always rewrite
     # tractable functions in terms of familiar intractable ones.
diff --git a/sympy/series/limits.py b/sympy/series/limits.py
index da6dd10501..faebd5da6b 100644
--- a/sympy/series/limits.py
+++ b/sympy/series/limits.py
@@ -3,7 +3,7 @@
 from sympy.core.exprtools import factor_terms
 from sympy.core.numbers import Float, _illegal
 from sympy.functions.combinatorial.factorials import factorial
-from sympy.functions.elementary.complexes import (Abs, sign, arg)
+from sympy.functions.elementary.complexes import (Abs, sign, arg, re)
 from sympy.functions.elementary.exponential import (exp, log)
 from sympy.functions.special.gamma_functions import gamma
 from sympy.polys import PolynomialError, factor
@@ -73,8 +73,8 @@ def heuristics(e, z, z0, dir):
     """
 
     rv = None
-    if abs(z0) is S.Infinity:
-        rv = limit(e.subs(z, 1/z), z, S.Zero, "+" if z0 is S.Infinity else "-")
+    if z0 is S.Infinity:
+        rv = limit(e.subs(z, 1/z), z, S.Zero, "+")
         if isinstance(rv, Limit):
             return
     elif e.is_Mul or e.is_Add or e.is_Pow or e.is_Function:
@@ -230,6 +230,13 @@ def doit(self, **hints):
             raise NotImplementedError("Limits at complex "
                                     "infinity are not implemented")
 
+        if z0.is_infinite:
+            cdir = sign(z0)
+            cdir = cdir/abs(cdir)
+            e = e.subs(z, cdir*z)
+            dir = "-"
+            z0 = S.Infinity
+
         if hints.get('deep', True):
             e = e.doit(**hints)
             z = z.doit(**hints)
@@ -289,7 +296,7 @@ def set_signs(expr):
 
 
         if e.is_meromorphic(z, z0):
-            if abs(z0) is S.Infinity:
+            if z0 is S.Infinity:
                 newe = e.subs(z, 1/z)
                 # cdir changes sign as oo- should become 0+
                 cdir = -cdir
@@ -311,7 +318,7 @@ def set_signs(expr):
                 else:
                     return S.ComplexInfinity
 
-        if abs(z0) is S.Infinity:
+        if z0 is S.Infinity:
             if e.is_Mul:
                 e = factor_terms(e)
             newe = e.subs(z, 1/z)
@@ -329,6 +336,12 @@ def set_signs(expr):
                 r = self.pow_heuristics(e)
                 if r is not None:
                     return r
+            try:
+                coeff = newe.as_leading_term(z, cdir=cdir)
+                if coeff != newe and coeff.has(exp):
+                    return gruntz(coeff, z, 0, "-" if re(cdir).is_negative else "+")
+            except (ValueError, NotImplementedError, PoleError):
+                pass
         else:
             if isinstance(coeff, AccumBounds) and ex == S.Zero:
                 return coeff
@@ -354,6 +367,8 @@ def set_signs(expr):
                             return S.Infinity*sign(coeff)*S.NegativeOne**ex
                         else:
                             return S.ComplexInfinity
+                else:
+                    raise NotImplementedError("Not sure of sign of %s" % ex)
 
         # gruntz fails on factorials but works with the gamma function
         # If no factorial term is present, e should remain unchanged.


Test patch:
diff --git a/sympy/concrete/tests/test_gosper.py b/sympy/concrete/tests/test_gosper.py
index fc81b70769..77b642a9b7 100644
--- a/sympy/concrete/tests/test_gosper.py
+++ b/sympy/concrete/tests/test_gosper.py
@@ -4,6 +4,7 @@
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol
 from sympy.functions.combinatorial.factorials import (binomial, factorial)
+from sympy.functions.elementary.exponential import (exp, log)
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.special.gamma_functions import gamma
 from sympy.polys.polytools import Poly
@@ -53,7 +54,7 @@ def test_gosper_sum():
     # issue 6033:
     assert gosper_sum(
         n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b)), \
-        (n, 0, m)).simplify() == -(a*b)**m*gamma(a + 1) \
+        (n, 0, m)).simplify() == -exp(m*log(a) + m*log(b))*gamma(a + 1) \
         *gamma(b + 1)/(gamma(a)*gamma(b)*gamma(a + m + 1)*gamma(b + m + 1)) \
         + 1/(gamma(a)*gamma(b))
 
diff --git a/sympy/functions/combinatorial/tests/test_comb_numbers.py b/sympy/functions/combinatorial/tests/test_comb_numbers.py
index e44a42da78..42afabb45a 100644
--- a/sympy/functions/combinatorial/tests/test_comb_numbers.py
+++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py
@@ -11,7 +11,7 @@
 from sympy.functions.elementary.complexes import (im, re)
 from sympy.functions.elementary.integers import floor
 from sympy.polys.polytools import cancel
-from sympy.series.limits import limit
+from sympy.series.limits import limit, Limit
 from sympy.series.order import O
 from sympy.functions import (
     bernoulli, harmonic, bell, fibonacci, tribonacci, lucas, euler, catalan,
@@ -335,6 +335,8 @@ def test_harmonic_evalf():
 def test_harmonic_rewrite():
     n = Symbol("n")
     m = Symbol("m", integer=True, positive=True)
+    x1 = Symbol("x1", positive=True)
+    x2 = Symbol("x2", negative=True)
 
     assert harmonic(n).rewrite(digamma) == polygamma(0, n + 1) + EulerGamma
     assert harmonic(n).rewrite(trigamma) == polygamma(0, n + 1) + EulerGamma
@@ -347,7 +349,9 @@ def test_harmonic_rewrite():
     assert expand_func(harmonic(n-4)) == harmonic(n) - 1/(n - 1) - 1/(n - 2) - 1/(n - 3) - 1/n
 
     assert harmonic(n, m).rewrite("tractable") == harmonic(n, m).rewrite(polygamma)
-    assert harmonic(n, x).rewrite("tractable") == zeta(x) - zeta(x, n+1)
+    assert harmonic(n, x1).rewrite("tractable") == harmonic(n, x1)
+    assert harmonic(n, x1 + 1).rewrite("tractable") == zeta(x1 + 1) - zeta(x1 + 1, n + 1)
+    assert harmonic(n, x2).rewrite("tractable") == zeta(x2) - zeta(x2, n + 1)
 
     _k = Dummy("k")
     assert harmonic(n).rewrite(Sum).dummy_eq(Sum(1/_k, (_k, 1, n)))
@@ -356,12 +360,18 @@ def test_harmonic_rewrite():
 
 def test_harmonic_calculus():
     y = Symbol("y", positive=True)
+    z = Symbol("z", negative=True)
     assert harmonic(x, 1).limit(x, 0) == 0
     assert harmonic(x, y).limit(x, 0) == 0
     assert harmonic(x, 1).series(x, y, 2) == \
             harmonic(y) + (x - y)*zeta(2, y + 1) + O((x - y)**2, (x, y))
     assert limit(harmonic(x, y), x, oo) == harmonic(oo, y)
-    assert limit(harmonic(x, y+1), x, oo) == zeta(y+1)
+    assert limit(harmonic(x, y + 1), x, oo) == zeta(y + 1)
+    assert limit(harmonic(x, y - 1), x, oo) == harmonic(oo, y - 1)
+    assert limit(harmonic(x, z), x, oo) == Limit(harmonic(x, z), x, oo, dir='-')
+    assert limit(harmonic(x, z + 1), x, oo) == oo
+    assert limit(harmonic(x, z + 2), x, oo) == harmonic(oo, z + 2)
+    assert limit(harmonic(x, z - 1), x, oo) == Limit(harmonic(x, z - 1), x, oo, dir='-')
 
 
 def test_euler():
diff --git a/sympy/functions/special/tests/test_bessel.py b/sympy/functions/special/tests/test_bessel.py
index b0d3089334..5effdefcbf 100644
--- a/sympy/functions/special/tests/test_bessel.py
+++ b/sympy/functions/special/tests/test_bessel.py
@@ -59,6 +59,18 @@ def test_bessely_leading_term():
     assert bessely(1, 2*sqrt(x)).as_leading_term(x) == -1/(pi*sqrt(x))
 
 
+def test_besseli_leading_term():
+    assert besseli(0, x).as_leading_term(x) == 1
+    assert besseli(1, sin(x)).as_leading_term(x) == x/2
+    assert besseli(1, 2*sqrt(x)).as_leading_term(x) == sqrt(x)
+
+
+def test_besselk_leading_term():
+    assert besselk(0, x).as_leading_term(x) == -log(x) - S.EulerGamma + log(2)
+    assert besselk(1, sin(x)).as_leading_term(x) == 1/x
+    assert besselk(1, 2*sqrt(x)).as_leading_term(x) == 1/(2*sqrt(x))
+
+
 def test_besselj_series():
     assert besselj(0, x).series(x) == 1 - x**2/4 + x**4/64 + O(x**6)
     assert besselj(0, x**(1.1)).series(x) == 1 + x**4.4/64 - x**2.2/4 + O(x**6)
@@ -82,6 +94,14 @@ def test_bessely_series():
     const = 2*S.EulerGamma/pi - 2*log(2)/pi + 2*log(x)/pi
     assert bessely(0, x).series(x, n=4) == const + x**2*(-log(x)/(2*pi)\
         + (2 - 2*S.EulerGamma)/(4*pi) + log(2)/(2*pi)) + O(x**4*log(x))
+    assert bessely(1, x).series(x, n=4) == -2/(pi*x) + x*(log(x)/pi - log(2)/pi - \
+        (1 - 2*S.EulerGamma)/(2*pi)) + x**3*(-log(x)/(8*pi) + \
+        (S(5)/2 - 2*S.EulerGamma)/(16*pi) + log(2)/(8*pi)) + O(x**4*log(x))
+    assert bessely(2, x).series(x, n=4) == -4/(pi*x**2) - 1/pi + x**2*(log(x)/(4*pi) - \
+        log(2)/(4*pi) - (S(3)/2 - 2*S.EulerGamma)/(8*pi)) + O(x**4*log(x))
+    assert bessely(3, x).series(x, n=4) == -16/(pi*x**3) - 2/(pi*x) - \
+        x/(4*pi) + x**3*(log(x)/(24*pi) - log(2)/(24*pi) - \
+        (S(11)/6 - 2*S.EulerGamma)/(48*pi)) + O(x**4*log(x))
     assert bessely(0, x**(1.1)).series(x, n=4) == 2*S.EulerGamma/pi\
         - 2*log(2)/pi + 2.2*log(x)/pi + x**2.2*(-0.55*log(x)/pi\
         + (2 - 2*S.EulerGamma)/(4*pi) + log(2)/(2*pi)) + O(x**4*log(x))
@@ -95,18 +115,59 @@ def test_bessely_series():
     assert bessely(0, log(1 + x)).series(x, n=3) == const\
         - x/pi + x**2*(-log(x)/(2*pi) + (2 - 2*S.EulerGamma)/(4*pi)\
         + log(2)/(2*pi) + 5/(12*pi)) + O(x**3*log(x))
-    assert bessely(1, sin(x)).series(x, n=4) == -(1/pi)*(1 - 2*S.EulerGamma)\
-        * (-x**3/12 + x/2) + x*(log(x)/pi - log(2)/pi) + x**3*(-7*log(x)\
-        / (24*pi) - 1/(6*pi) + (Rational(5, 2) - 2*S.EulerGamma)/(16*pi)\
-        + 7*log(2)/(24*pi)) + O(x**4*log(x))
-    assert bessely(1, 2*sqrt(x)).series(x, n=3) == sqrt(x)*(log(x)/pi \
-        - (1 - 2*S.EulerGamma)/pi) + x**Rational(3, 2)*(-log(x)/(2*pi)\
-        + (Rational(5, 2) - 2*S.EulerGamma)/(2*pi))\
-        + x**Rational(5, 2)*(log(x)/(12*pi)\
-        - (Rational(10, 3) - 2*S.EulerGamma)/(12*pi)) + O(x**3*log(x))
+    assert bessely(1, sin(x)).series(x, n=4) == -1/(pi*(-x**3/12 + x/2)) - \
+        (1 - 2*S.EulerGamma)*(-x**3/12 + x/2)/pi + x*(log(x)/pi - log(2)/pi) + \
+        x**3*(-7*log(x)/(24*pi) - 1/(6*pi) + (S(5)/2 - 2*S.EulerGamma)/(16*pi) +
+        7*log(2)/(24*pi)) + O(x**4*log(x))
+    assert bessely(1, 2*sqrt(x)).series(x, n=3) == -1/(pi*sqrt(x)) + \
+        sqrt(x)*(log(x)/pi - (1 - 2*S.EulerGamma)/pi) + x**(S(3)/2)*(-log(x)/(2*pi) + \
+        (S(5)/2 - 2*S.EulerGamma)/(2*pi)) + x**(S(5)/2)*(log(x)/(12*pi) - \
+        (S(10)/3 - 2*S.EulerGamma)/(12*pi)) + O(x**3*log(x))
     assert bessely(-2, sin(x)).series(x, n=4) == bessely(2, sin(x)).series(x, n=4)
 
 
+def test_besseli_series():
+    assert besseli(0, x).series(x) == 1 + x**2/4 + x**4/64 + O(x**6)
+    assert besseli(0, x**(1.1)).series(x) == 1 + x**4.4/64 + x**2.2/4 + O(x**6)
+    assert besseli(0, x**2 + x).series(x) == 1 + x**2/4 + x**3/2 + 17*x**4/64 + \
+        x**5/16 + O(x**6)
+    assert besseli(0, sqrt(x) + x).series(x, n=4) == 1 + x/4 + 17*x**2/64 + \
+        217*x**3/2304 + x**(S(3)/2)/2 + x**(S(5)/2)/16 + 25*x**(S(7)/2)/384 + O(x**4)
+    assert besseli(0, x/(1 - x)).series(x) == 1 + x**2/4 + x**3/2 + 49*x**4/64 + \
+        17*x**5/16 + O(x**6)
+    assert besseli(0, log(1 + x)).series(x) == 1 + x**2/4 - x**3/4 + 47*x**4/192 - \
+        23*x**5/96 + O(x**6)
+    assert besseli(1, sin(x)).series(x) == x/2 - x**3/48 - 47*x**5/1920 + O(x**6)
+    assert besseli(1, 2*sqrt(x)).series(x) == sqrt(x) + x**(S(3)/2)/2 + x**(S(5)/2)/12 + \
+        x**(S(7)/2)/144 + x**(S(9)/2)/2880 + x**(S(11)/2)/86400 + O(x**6)
+    assert besseli(-2, sin(x)).series(x, n=4) == besseli(2, sin(x)).series(x, n=4)
+
+
+def test_besselk_series():
+    const = log(2) - S.EulerGamma - log(x)
+    assert besselk(0, x).series(x, n=4) == const + \
+        x**2*(-log(x)/4 - S.EulerGamma/4 + log(2)/4 + S(1)/4) + O(x**4*log(x))
+    assert besselk(1, x).series(x, n=4) == 1/x + x*(log(x)/2 - log(2)/2 - \
+        S(1)/4 + S.EulerGamma/2) + x**3*(log(x)/16 - S(5)/64 - log(2)/16 + \
+        S.EulerGamma/16) + O(x**4*log(x))
+    assert besselk(2, x).series(x, n=4) == 2/x**2 - S(1)/2 + x**2*(-log(x)/8 - \
+        S.EulerGamma/8 + log(2)/8 + S(3)/32) + O(x**4*log(x))
+    assert besselk(0, x**(1.1)).series(x, n=4) == log(2) - S.EulerGamma - \
+        1.1*log(x) + x**2.2*(-0.275*log(x) - S.EulerGamma/4 + \
+        log(2)/4 + S(1)/4) + O(x**4*log(x))
+    assert besselk(0, x**2 + x).series(x, n=4) == const + \
+        (2 - 2*S.EulerGamma)*(x**3/4 + x**2/8) - x + x**2*(-log(x)/4 + \
+        log(2)/4 + S(1)/2) + x**3*(-log(x)/2 - S(7)/12 + log(2)/2) + O(x**4*log(x))
+    assert besselk(0, x/(1 - x)).series(x, n=3) == const - x + x**2*(-log(x)/4 - \
+        S(1)/4 - S.EulerGamma/4 + log(2)/4) + O(x**3*log(x))
+    assert besselk(0, log(1 + x)).series(x, n=3) == const + x/2 + \
+        x**2*(-log(x)/4 - S.EulerGamma/4 + S(1)/24 + log(2)/4) + O(x**3*log(x))
+    assert besselk(1, 2*sqrt(x)).series(x, n=3) == 1/(2*sqrt(x)) + \
+        sqrt(x)*(log(x)/2 - S(1)/2 + S.EulerGamma) + x**(S(3)/2)*(log(x)/4 - S(5)/8 + \
+        S.EulerGamma/2) + x**(S(5)/2)*(log(x)/24 - S(5)/36 + S.EulerGamma/12) + O(x**3*log(x))
+    assert besselk(-2, sin(x)).series(x, n=4) == besselk(2, sin(x)).series(x, n=4)
+
+
 def test_diff():
     assert besselj(n, z).diff(z) == besselj(n - 1, z)/2 - besselj(n + 1, z)/2
     assert bessely(n, z).diff(z) == bessely(n - 1, z)/2 - bessely(n + 1, z)/2
diff --git a/sympy/functions/special/tests/test_error_functions.py b/sympy/functions/special/tests/test_error_functions.py
index 3d8200148b..7d5afd116b 100644
--- a/sympy/functions/special/tests/test_error_functions.py
+++ b/sympy/functions/special/tests/test_error_functions.py
@@ -627,7 +627,8 @@ def test_si():
     t = Symbol('t', Dummy=True)
     assert Si(x).rewrite(sinc) == Integral(sinc(t), (t, 0, x))
 
-    assert limit(Shi(x), x, S.NegativeInfinity) == -I*pi/2
+    assert limit(Shi(x), x, S.Infinity) == S.Infinity
+    assert limit(Shi(x), x, S.NegativeInfinity) == S.NegativeInfinity
 
 
 def test_ci():
diff --git a/sympy/functions/special/tests/test_gamma_functions.py b/sympy/functions/special/tests/test_gamma_functions.py
index 0a200c534b..4be65bb470 100644
--- a/sympy/functions/special/tests/test_gamma_functions.py
+++ b/sympy/functions/special/tests/test_gamma_functions.py
@@ -618,6 +618,11 @@ def test_polygamma_expansion():
         2*x**3 + 3*x**4 + 2*x**5 - x**7 + 4*x**9/3 + O(x**11)
 
 
+def test_polygamma_leading_term():
+     expr = -log(1/x) + polygamma(0, 1 + 1/x) + S.EulerGamma
+     assert expr.as_leading_term(x, logx=-y) == S.EulerGamma
+
+
 def test_issue_8657():
     n = Symbol('n', negative=True, integer=True)
     m = Symbol('m', integer=True)
diff --git a/sympy/functions/special/tests/test_hyper.py b/sympy/functions/special/tests/test_hyper.py
index 716deac580..dc0eba0253 100644
--- a/sympy/functions/special/tests/test_hyper.py
+++ b/sympy/functions/special/tests/test_hyper.py
@@ -355,8 +355,6 @@ def test_limits():
     k, x = symbols('k, x')
     assert hyper((1,), (Rational(4, 3), Rational(5, 3)), k**2).series(k) == \
            1 + 9*k**2/20 + 81*k**4/1120 + O(k**6) # issue 6350
-    assert limit(meijerg((), (), (1,), (0,), -x), x, 0) == \
-            meijerg(((), ()), ((1,), (0,)), 0) # issue 6052
 
     # https://github.com/sympy/sympy/issues/11465
     assert limit(1/hyper((1, ), (1, ), x), x, 0) == 1
diff --git a/sympy/series/tests/test_gruntz.py b/sympy/series/tests/test_gruntz.py
index efcbca86c7..afb63b5977 100644
--- a/sympy/series/tests/test_gruntz.py
+++ b/sympy/series/tests/test_gruntz.py
@@ -5,6 +5,7 @@
 from sympy.functions.elementary.exponential import (exp, log)
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (acot, atan, cos, sin)
+from sympy.functions.elementary.complexes import sign as _sign
 from sympy.functions.special.error_functions import (Ei, erf)
 from sympy.functions.special.gamma_functions import (digamma, gamma, loggamma)
 from sympy.functions.special.zeta_functions import zeta
@@ -403,7 +404,7 @@ def test_gruntz_I():
     assert gruntz(I*x, x, oo) == I*oo
     assert gruntz(y*I*x, x, oo) == y*I*oo
     assert gruntz(y*3*I*x, x, oo) == y*I*oo
-    assert gruntz(y*3*sin(I)*x, x, oo) == y*I*oo
+    assert gruntz(y*3*sin(I)*x, x, oo).simplify().rewrite(_sign) == _sign(y)*I*oo
 
 
 def test_issue_4814():
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 8da96d0515..baf78e7e61 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -15,12 +15,14 @@
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.trigonometric import (acos, acot, acsc, asec, asin,
                                                       atan, cos, cot, csc, sec, sin, tan)
-from sympy.functions.special.bessel import (besselj, besselk)
+from sympy.functions.special.bessel import (besseli, bessely, besselj, besselk)
 from sympy.functions.special.error_functions import (Ei, erf, erfc, erfi, fresnelc, fresnels)
 from sympy.functions.special.gamma_functions import (digamma, gamma, uppergamma)
+from sympy.functions.special.hyper import meijerg
 from sympy.integrals.integrals import (Integral, integrate)
 from sympy.series.limits import (Limit, limit)
 from sympy.simplify.simplify import (logcombine, simplify)
+from sympy.simplify.hyperexpand import hyperexpand
 
 from sympy.calculus.accumulationbounds import AccumBounds
 from sympy.core.mul import Mul
@@ -346,6 +348,31 @@ def test_series_AccumBounds():
     assert limit(e, x, oo) == AccumBounds(0, oo)
 
 
+def test_bessel_functions_at_infinity():
+    # Pull Request 23844 implements limits for all bessel and modified bessel
+    # functions approaching infinity along any direction i.e. abs(z0) tends to oo
+
+    assert limit(besselj(1, x), x, oo) == 0
+    assert limit(besselj(1, x), x, -oo) == 0
+    assert limit(besselj(1, x), x, I*oo) == oo*I
+    assert limit(besselj(1, x), x, -I*oo) == -oo*I
+    assert limit(bessely(1, x), x, oo) == 0
+    assert limit(bessely(1, x), x, -oo) == 0
+    assert limit(bessely(1, x), x, I*oo) == -oo
+    assert limit(bessely(1, x), x, -I*oo) == -oo
+    assert limit(besseli(1, x), x, oo) == oo
+    assert limit(besseli(1, x), x, -oo) == -oo
+    assert limit(besseli(1, x), x, I*oo) == 0
+    assert limit(besseli(1, x), x, -I*oo) == 0
+    assert limit(besselk(1, x), x, oo) == 0
+    assert limit(besselk(1, x), x, -oo) == -oo*I
+    assert limit(besselk(1, x), x, I*oo) == 0
+    assert limit(besselk(1, x), x, -I*oo) == 0
+
+    # test issue 14874
+    assert limit(besselk(0, x), x, oo) == 0
+
+
 @XFAIL
 def test_doit2():
     f = Integral(2 * x, x)
@@ -461,6 +488,11 @@ def test_issue_3934():
 
 
 def test_calculate_series():
+    # NOTE
+    # The calculate_series method is being deprecated and is no longer responsible
+    # for result being returned. The mrv_leadterm function now uses simple leadterm
+    # calls rather than calculate_series.
+
     # needs gruntz calculate_series to go to n = 32
     assert limit(x**Rational(77, 3)/(1 + x**Rational(77, 3)), x, oo) == 1
     # needs gruntz calculate_series to go to n = 128
@@ -529,7 +561,7 @@ def test_factorial():
     # see Stirling's approximation:
     # https://en.wikipedia.org/wiki/Stirling's_approximation
     assert limit(f/(sqrt(2*pi*x)*(x/E)**x), x, oo) == 1
-    assert limit(f, x, -oo) == factorial(-oo)
+    assert limit(f, x, -oo) == gamma(-oo)
 
 
 def test_issue_6560():
@@ -630,6 +662,19 @@ def test_issue_4503():
         exp(x)/(2*sqrt(exp(x) + 1))
 
 
+def test_issue_6052():
+    G = meijerg((), (), (1,), (0,), -x)
+    g = hyperexpand(G)
+    assert limit(g, x, 0, '+-') == 0
+    assert limit(g, x, oo) == -oo
+
+
+def test_issue_7224():
+    expr = sqrt(x)*besseli(1,sqrt(8*x))
+    assert limit(x*diff(expr, x, x)/expr, x, 0) == 2
+    assert limit(x*diff(expr, x, x)/expr, x, 1).evalf(n=2) == 2.0
+
+
 def test_issue_8208():
     assert limit(n**(Rational(1, 1e9) - 1), n, oo) == 0
 
@@ -914,10 +959,6 @@ def test_issue_14811():
     assert limit(((1 + ((S(2)/3)**(x + 1)))**(2**x))/(2**((S(4)/3)**(x - 1))), x, oo) == oo
 
 
-def test_issue_14874():
-    assert limit(besselk(0, x), x, oo) == 0
-
-
 def test_issue_16222():
     assert limit(exp(x), x, 1000000000) == exp(1000000000)
 
@@ -1046,6 +1087,11 @@ def test_issue_16708():
     assert limit((B*ti*vi - sqrt(m)*sqrt(-2*B*d*vi + m*(vi)**2) + m*vi)/(B*vi), B, 0) == (d + ti*vi)/vi
 
 
+def test_issue_19154():
+    assert limit(besseli(1, 3 *x)/(x *besseli(1, x)**3), x , oo) == 2*sqrt(3)*pi/3
+    assert limit(besseli(1, 3 *x)/(x *besseli(1, x)**3), x , -oo) == -2*sqrt(3)*pi/3
+
+
 def test_issue_19453():
     beta = Symbol("beta", positive=True)
     h = Symbol("h", positive=True)
@@ -1206,6 +1252,14 @@ def test_issue_22220():
     assert limit(expr, x, pi, '+') == 2*I*pi
     assert limit(expr, x, pi, '-') == 0
 
+
+def test_issue_22334():
+    k, n  = symbols('k, n', positive=True)
+    assert limit((n+1)**k/((n+1)**(k+1) - (n)**(k+1)), n, oo) == 1/(k + 1)
+    assert limit((n+1)**k/((n+1)**(k+1) - (n)**(k+1)).expand(), n, oo) == 1/(k + 1)
+    assert limit((n+1)**k/(n*(-n**k + (n + 1)**k) + (n + 1)**k), n, oo) == 1/(k + 1)
+
+
 def test_sympyissue_22986():
     assert limit(acosh(1 + 1/x)*sqrt(x), x, oo) == sqrt(2)
 
diff --git a/sympy/series/tests/test_limitseq.py b/sympy/series/tests/test_limitseq.py
index 2cd0d4b6c6..362bb0397f 100644
--- a/sympy/series/tests/test_limitseq.py
+++ b/sympy/series/tests/test_limitseq.py
@@ -131,13 +131,19 @@ def test_issue_9308():
 
 def test_issue_10382():
     n = Symbol('n', integer=True)
-    assert limit_seq(fibonacci(n+1)/fibonacci(n), n) == S.GoldenRatio
+    assert limit_seq(fibonacci(n+1)/fibonacci(n), n).together() == S.GoldenRatio
 
 
 def test_issue_11672():
     assert limit_seq(Rational(-1, 2)**n, n) == 0
 
 
+def test_issue_14196():
+    k, n  = symbols('k, n', positive=True)
+    m = Symbol('m')
+    assert limit_seq(Sum(m**k, (m, 1, n)).doit()/(n**(k + 1)), n) == 1/(k + 1)
+
+
 def test_issue_16735():
     assert limit_seq(5**n/factorial(n), n) == 0
 
diff --git a/sympy/series/tests/test_series.py b/sympy/series/tests/test_series.py
index 7eca2b14d4..dd6eff3323 100644
--- a/sympy/series/tests/test_series.py
+++ b/sympy/series/tests/test_series.py
@@ -6,12 +6,12 @@
 from sympy.functions.elementary.exponential import (LambertW, exp, log)
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (atan, cos, sin)
-from sympy.integrals.integrals import Integral
+from sympy.functions.special.gamma_functions import gamma
+from sympy.integrals.integrals import Integral, integrate
 from sympy.series.order import O
 from sympy.series.series import series
 from sympy.abc import x, y, n, k
 from sympy.testing.pytest import raises
-from sympy.series.gruntz import calculate_series
 
 
 def test_sin():
@@ -88,6 +88,15 @@ def test_issue_5223():
     assert exp(sin(x)*log(x)).series(n=2) == 1 + x*log(x) + O(x**2*log(x)**2)
 
 
+def test_issue_6350():
+    expr = integrate(exp(k*(y**3 - 3*y)), (y, 0, oo), conds='none')
+    assert expr.series(k, 0, 3) == -(-1)**(S(2)/3)*sqrt(3)*gamma(S(1)/3)**2*gamma(S(2)/3)/(6*pi*k**(S(1)/3)) - \
+        sqrt(3)*k*gamma(-S(2)/3)*gamma(-S(1)/3)/(6*pi) - \
+        (-1)**(S(1)/3)*sqrt(3)*k**(S(1)/3)*gamma(-S(1)/3)*gamma(S(1)/3)*gamma(S(2)/3)/(6*pi) - \
+        (-1)**(S(2)/3)*sqrt(3)*k**(S(5)/3)*gamma(S(1)/3)**2*gamma(S(2)/3)/(4*pi) - \
+        (-1)**(S(1)/3)*sqrt(3)*k**(S(7)/3)*gamma(-S(1)/3)*gamma(S(1)/3)*gamma(S(2)/3)/(8*pi) + O(k**3)
+
+
 def test_issue_11313():
     assert Integral(cos(x), x).series(x) == sin(x).series(x)
     assert Derivative(sin(x), x).series(x, n=3).doit() == cos(x).series(x, n=3)
@@ -185,11 +194,6 @@ def test_x_is_base_detection():
     assert eq.series() == x**Rational(4, 3)
 
 
-def test_sin_power():
-    e = sin(x)**1.2
-    assert calculate_series(e, x) == x**1.2
-
-
 def test_issue_7203():
     assert series(cos(x), x, pi, 3) == \
         -1 + (x - pi)**2/2 + O((x - pi)**3, (x, pi))