Repository: sympy/sympy

Issue:
Bug in integration solver
For calculating the integration of pdf of cauchy distribution I tried sympy, but it produced wrong result.

`integrate(1/(pi*(1+(x-a)**2)),(x,-oo,oo))` 
the **expected** output was **1** but it **produced output 0 which is wrong**.

![image](https://user-images.githubusercontent.com/55033230/124940822-ef6a4e80-e027-11eb-8772-214cbb31fd8a.png)

![image](https://user-images.githubusercontent.com/55033230/124941424-67387900-e028-11eb-833e-e1b44d074e06.png)
![image](https://user-images.githubusercontent.com/55033230/124941649-95b65400-e028-11eb-99cf-26861bdab0f3.png)




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
  * Refactored leading term method for log Function. ([#23592](https://github.com/sympy/sympy/pull/23592) by [@anutosh491](https://github.com/anutosh491))

  * Added logx parameter to leading term methods of type add and power. ([#23592](https://github.com/sympy/sympy/pull/23592) by [@anutosh491](https://github.com/anutosh491))

  * Introduced arg_flag in Limit.doit.set_signs. ([#23592](https://github.com/sympy/sympy/pull/23592) by [@anutosh491](https://github.com/anutosh491))

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

<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
    This PR introduces a new algorithm/logic to compute leading terms of expressions of type `log` which has been broken down into 5 step explained below.
    1) In leading term methods the parameter `x` is expected to tend to `0` along the half-ray `x = cdir*t` where `t` is `real` and `positive`. Hence we define a positive dummy variable `t` and substitute `cdir*t` for `x` in the argument of `log: z = arg0.subs(x, cdir*t)`. 
    2) In the next step we calculate coefficient and exponent of the leading term of `z` through `c, e = z.leadterm(t)` to take decisions based on `c` and `e` for the incoming steps.
    3) In this step we handle the standard case of  `c, e == 1, 0`
    4) Next, we expand `log(c*t**e)` using `log(t) = log(x) - log(cdir)` or `logx - log(cdir)`.
    5) Finally we append `-2*I*pi` if `c` is `negative` and the imaginary part of `z` is `negative` for small `t > 0`

    #### 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
      * Refactored leading term method for log Function.
      * Added logx parameter to leading term methods of type add and power.
      * Introduced arg_flag in Limit.doit.set_signs.
    <!-- END RELEASE NOTES -->


</details><p>

anutosh491: Note : Limits/Series which need to be fixed . (These involve simply replacing the functions with their `log` rewrites and testing whether the new code is able to answer all of these).

1) Problems at `I*oo/ -I*oo` . Some handling is missing for these cases
```
>>> limit(atan(x), x, -I*oo)
-pi/2
>>> limit(atan(x).rewrite(log), x, -I*oo)     # expected pi/2
-pi/2
>>>
>>> limit(asinh(x), x, I*oo)
-oo
>>> limit(asinh(x).rewrite(log), x, I*oo)
log(zoo + oo*I)
>>> asinh(x).rewrite(log).series(x, I*oo)
log(zoo + oo*I)
>>>
>>> log(1 + x).series(x, x0=I*oo)    
O(1, (x, Limit(x + oo*I, x, 0)))    
>>> # expected answer for above
1/(5*x**5) - 1/(4*x**4) + 1/(3*x**3) - 1/(2*x**2) + 1/x + I*pi/2 + log(-I*x) + O(x**(-6), (x, oo*I)))
>>>
```
2) Limits/Series which are hang infinitely
```
>>> limit(acos(x).rewrite(log), x, 1)
>>> limit(acosh(x).rewrite(log), x, -oo)
>>> limit(asech(x).rewrite(log), x, -1)
>>> limit(asech(x).rewrite(log), x, 0, '-')
```
anutosh491: Also the second set of errors which hang infinitely do so because of `log._eval_nseries()` . For eg
```
>>> expr =  log(sqrt(1 - (x + 1)**2) + I*(x + 1))
>>> expr._eval_nseries(x, 0, None)
I*pi/2 + O(x**2)                # expected ans:  I*pi/2 + sqrt(2)*sqrt(x) + O(x**2)
```

jksuom: > im(z) != 0

It looks like this test does not suffice guarantee that not all terms of the lseries be real.
jksuom: The loop should be fixed by something like this:

        for i, term in enumerate(z.lseries(t)):
            if not term.is_real or i == 5:
                break
        if i < 5:
            coeff, _ = term.as_coeff_exponent(t)
            res += -2*I*S.Pi*Heaviside(-im(coeff), 0)

The test `im(z) != 0` can be improved but it is probably impossible to cover all cases, so limiting the loop count looks like the best solution. 
anutosh491: > > im(z) != 0
> 
> It looks like this test does not suffice to guarantee that not all terms of the lseries be real.

Yeah correct , wouldn't just using `im(z).is_negative` work here (I guess that's the main case we're trying to address here .. where we would have to subtract the `2*S.Pi*I` from `res`).  Not saying that's the permanent fix here but addresses the relevant problems here quickly 
```
>>> expr = Ei(-2*log(x))/2 - log(log(x))/(2*x**2)
>>> limit(expr , x, 1, '-')
EulerGamma/2 + log(2)/2 - I*pi/2
>>> limit(acosh(x).rewrite(log), x, -oo)
oo
>>> limit(asech(x).rewrite(log), x, -1)
I*pi
>>> limit(asech(x).rewrite(log), x, 0, '-')
oo
```
jksuom: > wouldn't just using im(z).is_negative work here

That does not cover the case where the sign of `im(z)` is unknown, for example `z = -1 + a*t` where `a` is a symbol.
anutosh491: Fixed a failing test 
```
>>> integrate(log(log(x)) / x**3, [x, 1, oo])
>>> # The above test from test_integrals.py would be calling the following
>>> limit(Ei(-2*log(x))/2 - log(log(x))/(2*x**2), x, oo)
NotImplementedError: Don't know how to calculate the mrv of 'O(_w, (_w, oo))'
```
And the workflow would involved leading term methods of `Ei` and `Chi` function , both of which were implemented wrongly ! For eg
```
>>> # On master
>>> limit(Ei(x), x, 0)
EulerGamma
>>> limit(Ei(x) - log(x), x, 0)
EulerGamma
>>> limit(Chi(x), x, 0)
EulerGamma
>>>
>>> # On branch
>>> limit(Ei(x), x, 0)
-oo
>>> limit(Ei(x) - log(x), x, 0)
EulerGamma
>>> limit(Chi(x), x, 0)
-oo
```
[Reference link](https://www.wolframalpha.com/input?i=limit+calculator&assumption=%7B%22F%22%2C+%22Limit%22%2C+%22limit%22%7D+-%3E%220%22&assumption=%7B%22F%22%2C+%22Limit%22%2C+%22limitfunction%22%7D+-%3E%22Ei%28x%29%22)

This also closes closely related issues (6682, 23596) for which I've added tests !
anutosh491: In the latest commit I've fixed leading term methods at points lying on branch cuts for all inverse trigonometric functions . I've added couple of TODO's which are remaining from our side
```
    # TODO : acos(x + 1).as_leading_term(x) hangs infinitely
    # assert acos(x + 1).as_leading_term(x) == sqrt(2)*sqrt(x)

    # TODO : atan(x - I).as_leading_term(x, cdir = -1) should be equal to
    # I*(log(x) + I*pi/2)/2 but currently gives I*(log(x) - I*pi/2)/2
    # assert atan(x - I).as_leading_term(x, cdir=-1) == I*(log(x) + I*pi/2)/2
```
Graph being 
![image](https://user-images.githubusercontent.com/87052487/173239375-44706a93-e7a8-47f3-ab78-b8d045e345b3.png)
So leading term from right side is fine but left side I've added as TODO 

Some basic limits which weren't responsive enough now work
```
>>> limit(atan(x), x, I)  # now I*oo
Limit(atan(x), x, I)
>>> limit(acot(x), x, I)  # now  -I*oo
Limit(acot(x), x, I)
>>> limit(asec(x), x, 0)   # now I*oo
zoo  
>>> limit(acsc(x), x, 0, '+')  # now -I*oo
zoo
```

anutosh491: Changes introduces in last few commits .
1) After fixing issues related to `add._eval_as_leading_term()`, few tests started failing in `test_bessel.py` . Good that they failed as leading term for `bessely` was buggy which has been fixed now . For eg
```
>>> # On master  All results here are wrong due to minute errors in leading term method
>>> bessely(0, x).as_leading_term(x)
(2*log(x) - 2*log(2))/pi
>>> limit(bessely(1, x), x, 0)
0
>>> limit(bessely(2, x), x, 0)
oo
>>> limit(bessely(3, x), x, 0)
oo
>>> limit(bessely(-1, x), x, 0)
0
>>> limit(bessely(-2, x), x, 0)
oo
>>> limit(bessely(-3, x), x, 0)
-oo
>>>
>>> # On branch
>>> bessely(0, x).as_leading_term(x)
(2*log(x) - 2*log(2) + 2*EulerGamma)/pi
>>> limit(bessely(1, x), x, 0)
-oo
>>> limit(bessely(2, x), x, 0)
-oo
>>> limit(bessely(3, x), x, 0)
-oo
>>> limit(bessely(-1, x), x, 0)
oo
>>> limit(bessely(-2, x), x, 0)
-oo
>>> limit(bessely(-3, x), x, 0)
oo
````
2) Another thing we might have to consider is that after introducing this change in the leading term method for add
```
        if logx is None:
            _logx = Dummy('logx')
            lgx = _logx
        else:
            lgx = logx
        leading_terms = [t.as_leading_term(x, logx=lgx, cdir=cdir) for t in expr.args]
```
Expressions involving multiple log arguments might be evaluated differently . we won't get wrong results as such but consider a test in `test_order.py`
```
>>> O(log(x) + 1/log(x))
O(log(x))
>>> O(log(x) + 1/log(x))
O((log(x)**2 + 1)/log(x))
```
The second result isn't wrong either and is the only test concerning the change as of now , we could maybe look into this more carefully if required !

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
       before           after         ratio
     [77f1d79c]       [d0134b30]
     <sympy-1.10.1^0>                 
+         126±4ms          235±3ms     1.87  sum.TimeSum.time_doit

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
anutosh491: I feel we should be able to answer #21721 through the changes introduced here, we are somewhat close . Boils down to this step 

```
>>> (-I*log(-a - I + 1/x)+ I*log(-a + I + 1/x))
-I*log(-a - I + 1/x) + I*log(-a + I + 1/x)
>>> _.as_leading_term(x, cdir =-1)   # expected 0
2*pi*Heaviside(im(a) - 1, 0) - 2*pi*Heaviside(im(a) + 1, 0)
>>>
>>> expr1 = -I*log(-a - I + 1/x)
>>> expr2 = I*log(-a + I + 1/x)
>>> expr1.as_leading_term(x, cdir=-1).expand()
I*log(x) - 2*pi*Heaviside(im(a) + 1, 0) + 2*pi
>>> expr2.as_leading_term(x, cdir=-1).expand()
-I*log(x) + 2*pi*Heaviside(im(a) - 1, 0) - 2*pi
```
The expected result here is `0` , the limit/leading_term doesn't really depend on `a`. I guess the series shows that
```
>>> expr = (-I*log(-a - I + 1/x)+ I*log(-a + I + 1/x))
>>> expr.series(x, 0, cdir=-1)
-2*x + x**3*(2/3 - 2*a**2) + x**4*(-2*a**3 + 2*a) + x**5*(-2*a**4 + 4*a**2 - 2/5) - 2*a*x**2 + O(x**6)
>>> _.as_leading_term(x, cdir =-1)
-2*x
```
Though the algorithm doesn't seem to make any obvious mistakes here, need to think how we can handle such cases !

Patch:
diff --git a/sympy/concrete/summations.py b/sympy/concrete/summations.py
index 830785990c..4890ebc3b5 100644
--- a/sympy/concrete/summations.py
+++ b/sympy/concrete/summations.py
@@ -505,7 +505,8 @@ def is_convergent(self):
 
         ### ------------- comparison test ------------- ###
         # 1/(n**p*log(n)**q*log(log(n))**r) comparison
-        n_log_test = order.expr.match(1/(sym**p*log(sym)**q*log(log(sym))**r))
+        n_log_test = (order.expr.match(1/(sym**p*log(1/sym)**q*log(-log(1/sym))**r)) or
+                      order.expr.match(1/(sym**p*(-log(1/sym))**q*log(-log(1/sym))**r)))
         if n_log_test is not None:
             if (n_log_test[p] > 1 or
                 (n_log_test[p] == 1 and n_log_test[q] > 1) or
diff --git a/sympy/core/add.py b/sympy/core/add.py
index 73c5a18f55..d3223dc316 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -993,12 +993,16 @@ 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.series.order import Order
         from sympy.functions.elementary.exponential import log
         from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
         from .function import expand_mul
 
-        old = self
+        o = self.getO()
+        if o is None:
+            o = Order(0)
+        old = self.removeO()
 
         if old.has(Piecewise):
             old = piecewise_fold(old)
@@ -1017,7 +1021,8 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
 
         infinite = [t for t in expr.args if t.is_infinite]
 
-        leading_terms = [t.as_leading_term(x, logx=logx, cdir=cdir) for t in expr.args]
+        _logx = Dummy('logx') if logx is None else logx
+        leading_terms = [t.as_leading_term(x, logx=_logx, cdir=cdir) for t in expr.args]
 
         min, new_expr = Order(0), 0
 
@@ -1033,6 +1038,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         except TypeError:
             return expr
 
+        if logx is None:
+            new_expr = new_expr.subs(_logx, log(x))
+
         is_zero = new_expr.is_zero
         if is_zero is None:
             new_expr = new_expr.trigsimp().cancel()
@@ -1040,16 +1048,19 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if is_zero is True:
             # simple leading term analysis gave us cancelled terms but we have to send
             # back a term, so compute the leading term (via series)
-            n0 = min.getn()
+            try:
+                n0 = min.getn()
+            except NotImplementedError:
+                n0 = 1
             res = Order(1)
             incr = S.One
             while res.is_Order:
-                res = old._eval_nseries(x, n=n0+incr, logx=None, cdir=cdir).cancel().powsimp().trigsimp()
+                res = old._eval_nseries(x, n=n0+incr, logx=logx, cdir=cdir).cancel().powsimp().trigsimp()
                 incr *= 2
             return res.as_leading_term(x, logx=logx, cdir=cdir)
 
         elif new_expr is S.NaN:
-            return old.func._from_args(infinite)
+            return old.func._from_args(infinite) + o
 
         else:
             return new_expr
diff --git a/sympy/core/mul.py b/sympy/core/mul.py
index 3c010d8b85..bac30d0b37 100644
--- a/sympy/core/mul.py
+++ b/sympy/core/mul.py
@@ -1902,7 +1902,7 @@ def coeff_exp(term, x):
 
         try:
             for t in self.args:
-                coeff, exp = t.leadterm(x, logx=logx)
+                coeff, exp = t.leadterm(x)
                 if not coeff.has(x):
                     ords.append((t, exp))
                 else:
@@ -1965,6 +1965,10 @@ def max_degree(e, x):
                 return res
 
         if res != self:
+            if (self - res).subs(x, 0) == S.Zero and n > 0:
+                lt = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
+                if lt == S.Zero:
+                    return res
             res += Order(x**n, x)
         return res
 
diff --git a/sympy/core/power.py b/sympy/core/power.py
index 018ae5d307..b34ae87f50 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1723,7 +1723,7 @@ def mul(d1, d2):
             return res
 
         try:
-            _, d = g.leadterm(x)
+            _, d = g.leadterm(x, logx=logx)
         except (ValueError, NotImplementedError):
             if limit(g/x**maxpow, x, 0) == 0:
                 # g has higher order zero
@@ -1734,10 +1734,10 @@ def mul(d1, d2):
             g = g.simplify()
             if g.is_zero:
                 return f**e
-            _, d = g.leadterm(x)
+            _, d = g.leadterm(x, logx=logx)
             if not d.is_positive:
                 g = ((b - f)/f).expand()
-                _, d = g.leadterm(x)
+                _, d = g.leadterm(x, logx=logx)
                 if not d.is_positive:
                     raise NotImplementedError()
 
@@ -1797,7 +1797,10 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             return lt.as_leading_term(x, logx=logx, cdir=cdir)
         else:
             from sympy.functions.elementary.complexes import im
-            f = b.as_leading_term(x, logx=logx, cdir=cdir)
+            try:
+                f = b.as_leading_term(x, logx=logx, cdir=cdir)
+            except PoleError:
+                return self
             if (not e.is_integer and f.is_constant() and f.is_real
                 and f.is_negative and im((b - f).dir(x, cdir)).is_negative):
                 return self.func(f, e) * exp(-2 * e * S.Pi * S.ImaginaryUnit)
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index a3994b3a04..0a36384748 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -507,6 +507,8 @@ def _eval_nseries(self, x, n, logx, cdir=0):
             nterms = ceiling(n/cf)
         exp_series = exp(t)._taylor(t, nterms)
         r = exp(arg0)*exp_series.subs(t, arg_series - arg0)
+        if r.subs(logx, log(x)) == self:
+            return r
         if cf and cf > 1:
             r += Order((arg_series - arg0)**n, x)/x**((cf-1)*n)
         else:
@@ -971,6 +973,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):
         #      IMPORTANT.
         from sympy.series.order import Order
         from sympy.simplify.simplify import logcombine
+        _logx = logx
         if not logx:
             logx = log(x)
         if self.args[0] == x:
@@ -1007,7 +1010,13 @@ def coeff_exp(term, x):
             while s.is_Order:
                 n += 1
                 s = arg.nseries(x, n=n, logx=logx)
-        a, b = s.removeO().leadterm(x)
+        if _logx and logx.has(x):
+            a, b = arg.as_leading_term(x, logx=logx), S.Zero
+        else:
+            try:
+                a, b = s.removeO().leadterm(x)
+            except (ValueError, NotImplementedError, PoleError):
+                a, b = s.removeO().as_leading_term(x), S.Zero
         p = cancel(s/(a*x**b) - 1).expand().powsimp()
         if p.has(exp):
             p = logcombine(p)
@@ -1015,7 +1024,20 @@ def coeff_exp(term, x):
             n = p.getn()
         _, d = coeff_exp(p, x)
         if not d.is_positive:
-            return log(a) + b*logx + Order(x**n, x)
+            res = log(a) + b*logx
+            _res = res
+            logflags = dict(deep=True, log=True, mul=False, power_exp=False,
+                power_base=False, multinomial=False, basic=False, force=True,
+                factor=False)
+            expr = self.expand(**logflags)
+            if (not a.could_extract_minus_sign() and
+                logx.could_extract_minus_sign()):
+                _res = _res.subs(-logx, -log(x)).expand(**logflags)
+            else:
+                _res = _res.subs(logx, log(x)).expand(**logflags)
+            if _res == expr:
+                return res
+            return res + Order(x**n, x)
 
         def mul(d1, d2):
             res = {}
@@ -1053,24 +1075,48 @@ def mul(d1, d2):
         return res + Order(x**n, x)
 
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
-        from sympy.calculus.util import AccumBounds
+        # NOTE
+        # Refer https://github.com/sympy/sympy/pull/23592 for more information
+        # on each of the following steps involved in this method.
         arg0 = self.args[0].together()
 
-        arg = arg0.as_leading_term(x, cdir=cdir)
-        x0 = arg0.subs(x, 0)
-        if (x0 is S.NaN and logx is None):
-            x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
-        if x0 in (S.NegativeInfinity, S.Infinity):
-            raise PoleError("Cannot expand %s around 0" % (self))
-        if x0 == 1:
-            return (arg0 - S.One).as_leading_term(x)
-        if cdir != 0:
-            cdir = arg0.dir(x, cdir)
-        if x0.is_real and x0.is_negative and im(cdir).is_negative:
-            return self.func(x0) - 2*I*S.Pi
-        if isinstance(arg, AccumBounds):
+        # STEP 1
+        t = Dummy('t', real=True, positive=True)
+        if cdir == 0:
+            cdir = 1
+        z = arg0.subs(x, cdir*t)
+
+        # STEP 2
+        try:
+            c, e = z.leadterm(t, logx=logx, cdir=1)
+        except ValueError:
+            arg = arg0.as_leading_term(x, logx=logx, cdir=cdir)
             return log(arg)
-        return self.func(arg)
+        if c.has(t):
+            c = c.subs(t, x/cdir)
+            if e != 0:
+                raise PoleError("Cannot expand %s around 0" % (self))
+            return log(c)
+
+        # STEP 3
+        if c == S.One and e == S.Zero:
+            return (arg0 - S.One).as_leading_term(x, logx=logx)
+
+        # STEP 4
+        res = log(c) - e*log(cdir)
+        logx = log(x) if logx is None else logx
+        res += e*logx
+
+        # STEP 5
+        if c.is_negative and im(z) != 0:
+            from sympy.functions.special.delta_functions import Heaviside
+            for i, term in enumerate(z.lseries(t)):
+                if not term.is_real or i == 5:
+                    break
+            if i < 5:
+                coeff, _ = term.as_coeff_exponent(t)
+                res += -2*I*S.Pi*Heaviside(-im(coeff), 0)
+        return res
 
 
 class LambertW(Function):
diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py
index dd237ba0b6..3b17ee8c8f 100644
--- a/sympy/functions/elementary/trigonometric.py
+++ b/sympy/functions/elementary/trigonometric.py
@@ -2318,8 +2318,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         x0 = arg.subs(x, 0).cancel()
         if x0.is_zero:
             return arg.as_leading_term(x)
-        if x0 is S.ComplexInfinity:
-            return S.ImaginaryUnit*log(arg.as_leading_term(x))
+        # Handling Branch cuts (-oo, -1) U (1, oo)
+        if x0 in (-S.One, S.One, S.ComplexInfinity):
+            return self.rewrite(log)._eval_as_leading_term(x, logx=logx, cdir=cdir).expand()
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 < S.NegativeOne:
@@ -2374,6 +2375,7 @@ def _eval_rewrite_as_atan(self, x, **kwargs):
 
     def _eval_rewrite_as_log(self, x, **kwargs):
         return -S.ImaginaryUnit*log(S.ImaginaryUnit*x + sqrt(1 - x**2))
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
 
     def _eval_rewrite_as_acot(self, arg, **kwargs):
         return 2*acot((1 + sqrt(1 - arg**2))/arg)
@@ -2518,10 +2520,11 @@ def taylor_term(n, x, *previous_terms):
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
+        # Handling Branch cuts (-oo, -1) U (1, oo)
         if x0 == 1:
             return sqrt(2)*sqrt((S.One - arg).as_leading_term(x))
-        if x0 is S.ComplexInfinity:
-            return S.ImaginaryUnit*log(arg.as_leading_term(x))
+        if x0 in (-S.One, S.ComplexInfinity):
+            return self.rewrite(log)._eval_as_leading_term(x, logx=logx, cdir=cdir)
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 < S.NegativeOne:
@@ -2578,6 +2581,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # acos
     def _eval_rewrite_as_log(self, x, **kwargs):
         return pi/2 + S.ImaginaryUnit*\
             log(S.ImaginaryUnit*x + sqrt(1 - x**2))
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
 
     def _eval_rewrite_as_asin(self, x, **kwargs):
         return pi/2 - asin(x)
@@ -2745,8 +2749,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         x0 = arg.subs(x, 0).cancel()
         if x0.is_zero:
             return arg.as_leading_term(x)
-        if x0 is S.ComplexInfinity:
-            return acot(1/arg)._eval_as_leading_term(x, cdir=cdir)
+        # Handling Branch cuts (-I*oo, -I) U (I, I*oo)
+        if x0 in (-S.ImaginaryUnit, S.ImaginaryUnit, S.ComplexInfinity):
+            return self.rewrite(log)._eval_as_leading_term(x, logx=logx, cdir=cdir).expand()
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if re(cdir) < 0 and re(x0).is_zero and im(x0) > S.One:
@@ -2773,6 +2778,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # atan
     def _eval_rewrite_as_log(self, x, **kwargs):
         return S.ImaginaryUnit/2*(log(S.One - S.ImaginaryUnit*x)
             - log(S.One + S.ImaginaryUnit*x))
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
 
     def _eval_aseries(self, n, args0, x, logx):
         if args0[0] is S.Infinity:
@@ -2944,10 +2950,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             return (1/arg).as_leading_term(x)
         if cdir != 0:
             cdir = arg.dir(x, cdir)
-        if x0.is_zero:
-            if re(cdir) < 0:
-                return self.func(x0) - pi
-            return self.func(x0)
+        # Handling Branch cuts [-I, I]
+        if x0 in (-S.ImaginaryUnit, S.ImaginaryUnit, S.Zero):
+            return self.rewrite(log)._eval_as_leading_term(x, logx=logx, cdir=cdir).expand()
         if re(cdir) > 0 and re(x0).is_zero and im(x0) > S.Zero and im(x0) < S.One:
             return self.func(x0) + pi
         if re(cdir) < 0 and re(x0).is_zero and im(x0) < S.Zero and im(x0) > S.NegativeOne:
@@ -2982,6 +2987,7 @@ def _eval_aseries(self, n, args0, x, logx):
     def _eval_rewrite_as_log(self, x, **kwargs):
         return S.ImaginaryUnit/2*(log(1 - S.ImaginaryUnit/x)
             - log(1 + S.ImaginaryUnit/x))
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
 
     def inverse(self, argindex=1):
         """
@@ -3116,10 +3122,11 @@ def inverse(self, argindex=1):
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
+        # Handling Branch cuts (-1, 1)
         if x0 == 1:
             return sqrt(2)*sqrt((arg - S.One).as_leading_term(x))
-        if x0.is_zero:
-            return S.ImaginaryUnit*log(arg.as_leading_term(x))
+        if x0 in (-S.One, S.Zero):
+            return self.rewrite(log)._eval_as_leading_term(x, logx=logx, cdir=cdir)
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 > S.Zero and x0 < S.One:
@@ -3170,6 +3177,7 @@ def _eval_is_extended_real(self):
 
     def _eval_rewrite_as_log(self, arg, **kwargs):
         return pi/2 + S.ImaginaryUnit*log(S.ImaginaryUnit/arg + sqrt(1 - 1/arg**2))
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
 
     def _eval_rewrite_as_asin(self, arg, **kwargs):
         return pi/2 - asin(1/arg)
@@ -3287,10 +3295,11 @@ def inverse(self, argindex=1):
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        if x0.is_zero:
-            return S.ImaginaryUnit*log(arg.as_leading_term(x))
+        # Handling Branch cuts (-1, 1)
+        if x0 in (-S.One, S.One, S.Zero):
+            return self.rewrite(log)._eval_as_leading_term(x, logx=logx, cdir=cdir).expand()
         if x0 is S.ComplexInfinity:
-            return arg.as_leading_term(x)
+            return (1/arg).as_leading_term(x)
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 > S.Zero and x0 < S.One:
@@ -3335,6 +3344,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # acsc
 
     def _eval_rewrite_as_log(self, arg, **kwargs):
         return -S.ImaginaryUnit*log(S.ImaginaryUnit/arg + sqrt(1 - 1/arg**2))
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
 
     def _eval_rewrite_as_asin(self, arg, **kwargs):
         return asin(1/arg)
diff --git a/sympy/functions/special/bessel.py b/sympy/functions/special/bessel.py
index 3d0df1ef80..3c2a07a76c 100644
--- a/sympy/functions/special/bessel.py
+++ b/sympy/functions/special/bessel.py
@@ -339,10 +339,10 @@ 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 - 1).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)
-        if x in arg.free_symbols:
+        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):
             return arg
         else:
             return self.func(nu, z.subs(x, 0).cancel())
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 832adc2175..b3b4572ff8 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -1212,10 +1212,15 @@ def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):
         return exp(z) * _eis(z)
 
     def _eval_as_leading_term(self, x, logx=None, cdir=0):
+        from sympy import re
         x0 = self.args[0].limit(x, 0)
+        arg = self.args[0].as_leading_term(x, cdir=cdir)
+        cdir = arg.dir(x, cdir)
         if x0.is_zero:
-            f = self._eval_rewrite_as_Si(*self.args)
-            return f._eval_as_leading_term(x, logx=logx, cdir=cdir)
+            c, e = arg.as_coeff_exponent(x)
+            logx = log(x) if logx is None else logx
+            return log(c) + e*logx + S.EulerGamma - (
+                S.ImaginaryUnit*pi if re(cdir).is_negative else S.Zero)
         return super()._eval_as_leading_term(x, logx=logx, cdir=cdir)
 
     def _eval_nseries(self, x, n, logx, cdir=0):
@@ -2251,7 +2256,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:
diff --git a/sympy/series/gruntz.py b/sympy/series/gruntz.py
index 4c2579c6d9..0e19e56b88 100644
--- a/sympy/series/gruntz.py
+++ b/sympy/series/gruntz.py
@@ -695,9 +695,9 @@ def gruntz(e, z, z0, dir="+"):
 
     # convert all limits to the limit z->oo; sign of z is handled in limitinf
     r = None
-    if z0 == oo:
+    if z0 in (oo, I*oo):
         e0 = e
-    elif z0 == -oo:
+    elif z0 in (-oo, -I*oo):
         e0 = e.subs(z, -z)
     else:
         if str(dir) == "-":
diff --git a/sympy/series/limits.py b/sympy/series/limits.py
index 69060db76e..1bac908dcd 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)
+from sympy.functions.elementary.complexes import (Abs, sign, arg)
 from sympy.functions.elementary.exponential import (exp, log)
 from sympy.functions.special.gamma_functions import gamma
 from sympy.polys import PolynomialError, factor
@@ -147,9 +147,9 @@ def __new__(cls, e, z, z0, dir="+"):
         z = sympify(z)
         z0 = sympify(z0)
 
-        if z0 is S.Infinity:
+        if z0 in (S.Infinity, S.ImaginaryUnit*S.Infinity):
             dir = "-"
-        elif z0 is S.NegativeInfinity:
+        elif z0 in (S.NegativeInfinity, S.ImaginaryUnit*S.NegativeInfinity):
             dir = "+"
 
         if(z0.has(z)):
@@ -249,16 +249,19 @@ def set_signs(expr):
             if newargs != expr.args:
                 expr = expr.func(*newargs)
             abs_flag = isinstance(expr, Abs)
+            arg_flag = isinstance(expr, arg)
             sign_flag = isinstance(expr, sign)
-            if abs_flag or sign_flag:
+            if abs_flag or sign_flag or arg_flag:
                 sig = limit(expr.args[0], z, z0, dir)
                 if sig.is_zero:
                     sig = limit(1/expr.args[0], z, z0, dir)
                 if sig.is_extended_real:
                     if (sig < 0) == True:
-                        return -expr.args[0] if abs_flag else S.NegativeOne
+                        return (-expr.args[0] if abs_flag else
+                                S.NegativeOne if sign_flag else S.Pi)
                     elif (sig > 0) == True:
-                        return expr.args[0] if abs_flag else S.One
+                        return (expr.args[0] if abs_flag else
+                                S.One if sign_flag else S.Zero)
             return expr
 
         if e.has(Float):


Test patch:
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index 9b53834298..e7e5717753 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -484,6 +484,27 @@ def test_log_apply_evalf():
     assert value.epsilon_eq(Float("0.58496250072115618145373"))
 
 
+def test_log_leading_term():
+    p = Symbol('p')
+
+    # Test for STEP 3
+    assert log(1 + x + x**2).as_leading_term(x, cdir=1) == x
+    # Test for STEP 4
+    assert log(2*x).as_leading_term(x, cdir=1) == log(x) + log(2)
+    assert log(2*x).as_leading_term(x, cdir=-1) == log(x) + log(2)
+    assert log(-2*x).as_leading_term(x, cdir=1, logx=p) == p + log(2) + I*pi
+    assert log(-2*x).as_leading_term(x, cdir=-1, logx=p) == p + log(2) - I*pi
+    # Test for STEP 5
+    assert log(-2*x + (3 - I)*x**2).as_leading_term(x, cdir=1) == log(x) + log(2) - I*pi
+    assert log(-2*x + (3 - I)*x**2).as_leading_term(x, cdir=-1) == log(x) + log(2) - I*pi
+    assert log(2*x + (3 - I)*x**2).as_leading_term(x, cdir=1) == log(x) + log(2)
+    assert log(2*x + (3 - I)*x**2).as_leading_term(x, cdir=-1) == log(x) + log(2) - 2*I*pi
+    assert log(-1 + x - I*x**2 + I*x**3).as_leading_term(x, cdir=1) == -I*pi
+    assert log(-1 + x - I*x**2 + I*x**3).as_leading_term(x, cdir=-1) == -I*pi
+    assert log(-1/(1 - x)).as_leading_term(x, cdir=1) == I*pi
+    assert log(-1/(1 - x)).as_leading_term(x, cdir=-1) == I*pi
+
+
 def test_log_nseries():
     assert log(x - 1)._eval_nseries(x, 4, None, I) == I*pi - x - x**2/2 - x**3/3 + O(x**4)
     assert log(x - 1)._eval_nseries(x, 4, None, -I) == -I*pi - x - x**2/2 - x**3/3 + O(x**4)
diff --git a/sympy/functions/elementary/tests/test_trigonometric.py b/sympy/functions/elementary/tests/test_trigonometric.py
index 2c77d50af4..bebbaa7fb0 100644
--- a/sympy/functions/elementary/tests/test_trigonometric.py
+++ b/sympy/functions/elementary/tests/test_trigonometric.py
@@ -900,7 +900,6 @@ def test_asin():
                 assert asin(sin(n*pi/d)) == n*pi/d
 
     assert asin(x).diff(x) == 1/sqrt(1 - x**2)
-    assert asin(1/x).as_leading_term(x) == I*log(1/x)
 
     assert asin(0.2, evaluate=False).is_real is True
     assert asin(-2).is_real is False
@@ -923,6 +922,13 @@ def test_asin_series():
     assert t5 == 3*x**5/40
     assert asin(x).taylor_term(7, x, t5, 0) == 5*x**7/112
 
+def test_asin_leading_term():
+    assert asin(x).as_leading_term(x) == x
+    assert asin(x + 1).as_leading_term(x) == pi/2
+    assert asin(x - 1).as_leading_term(x) == -pi/2
+    assert asin(1/x).as_leading_term(x, cdir=1) == I*log(x) + pi/2 - I*log(2)
+    assert asin(1/x).as_leading_term(x, cdir=-1) == -I*log(x) - 3*pi/2 + I*log(2)
+
 
 def test_asin_rewrite():
     assert asin(x).rewrite(log) == -I*log(I*x + sqrt(1 - x**2))
@@ -964,7 +970,6 @@ def test_acos():
     assert acos(2*I) == pi/2 - asin(2*I)
 
     assert acos(x).diff(x) == -1/sqrt(1 - x**2)
-    assert acos(1/x).as_leading_term(x) == I*log(1/x)
 
     assert acos(0.2).is_real is True
     assert acos(-2).is_real is False
@@ -983,6 +988,14 @@ def test_acos():
     assert acos(z).conjugate() != acos(conjugate(z))
 
 
+def test_acos_leading_term():
+    assert acos(x).as_leading_term(x) == pi/2
+    assert acos(x + 1).as_leading_term(x) == sqrt(2)*sqrt(-x)
+    assert acos(x - 1).as_leading_term(x) == pi
+    assert acos(1/x).as_leading_term(x, cdir=1) == -I*log(x) + I*log(2)
+    assert acos(1/x).as_leading_term(x, cdir=-1) == I*log(x) + 2*pi - I*log(2)
+
+
 def test_acos_series():
     assert acos(x).series(x, 0, 8) == \
         pi/2 - x - x**3/6 - 3*x**5/40 - 5*x**7/112 + O(x**8)
@@ -1038,7 +1051,6 @@ def test_atan():
 
     assert atan(oo) == pi/2
     assert atan(x).diff(x) == 1/(1 + x**2)
-    assert atan(1/x).as_leading_term(x) == pi/2
 
     assert atan(r).is_real is True
 
@@ -1084,6 +1096,16 @@ def test_atan_fdiff():
     raises(ArgumentIndexError, lambda: atan(x).fdiff(2))
 
 
+def test_atan_leading_term():
+    assert atan(x).as_leading_term(x) == x
+    assert atan(x + I).as_leading_term(x, cdir=1) == -I*log(x)/2 + pi/4 + I*log(2)/2
+    assert atan(x + I).as_leading_term(x, cdir=-1) == -I*log(x)/2 - 3*pi/4 + I*log(2)/2
+    assert atan(x - I).as_leading_term(x, cdir=1) == I*log(x)/2 + pi/4 - I*log(2)/2
+    assert atan(x - I).as_leading_term(x, cdir=-1) == I*log(x)/2 + pi/4 - I*log(2)/2
+    assert atan(1/x).as_leading_term(x, cdir=1) == pi/2
+    assert atan(1/x).as_leading_term(x, cdir=-1) == -pi/2
+
+
 def test_atan2():
     assert atan2.nargs == FiniteSet(2)
     assert atan2(0, 0) is S.NaN
@@ -1181,7 +1203,6 @@ def test_acot():
     assert acot(1/sqrt(3)) == pi/3
     assert acot(-1/sqrt(3)) == -pi/3
     assert acot(x).diff(x) == -1/(1 + x**2)
-    assert acot(1/x).as_leading_term(x) == x
 
     assert acot(r).is_extended_real is True
 
@@ -1214,6 +1235,15 @@ def test_acot_fdiff():
     assert acot(x).fdiff() == -1/(x**2 + 1)
     raises(ArgumentIndexError, lambda: acot(x).fdiff(2))
 
+def test_acot_leading_term():
+    assert acot(x).as_leading_term(x, cdir=1) == pi/2
+    assert acot(x).as_leading_term(x, cdir=-1) == -pi/2
+    assert acot(x + I).as_leading_term(x, cdir=1) == I*log(x)/2 + pi/4 - I*log(2)/2
+    assert acot(x + I).as_leading_term(x, cdir=-1) == I*log(x)/2 + pi/4 - I*log(2)/2
+    assert acot(x - I).as_leading_term(x, cdir=1) == -I*log(x)/2 + pi/4 + I*log(2)/2
+    assert acot(x - I).as_leading_term(x, cdir=-1) == -I*log(x)/2 - 3*pi/4 + I*log(2)/2
+    assert acot(1/x).as_leading_term(x) == x
+
 
 def test_attributes():
     assert sin(x).args == (x,)
@@ -1277,10 +1307,6 @@ def test_as_leading_term_issue_5272():
     assert cos(x).as_leading_term(x) == 1
     assert tan(x).as_leading_term(x) == x
     assert cot(x).as_leading_term(x) == 1/x
-    assert asin(x).as_leading_term(x) == x
-    assert acos(x).as_leading_term(x) == pi/2
-    assert atan(x).as_leading_term(x) == x
-    assert acot(x).as_leading_term(x) == pi/2
 
 
 def test_leading_terms():
@@ -1734,7 +1760,6 @@ def test_asec():
     assert asec(sqrt(2) - sqrt(6)) == pi*Rational(11, 12)
 
     assert asec(x).diff(x) == 1/(x**2*sqrt(1 - 1/x**2))
-    assert asec(x).as_leading_term(x) == I*log(x)
 
     assert asec(x).rewrite(log) == I*log(sqrt(1 - 1/x**2) + I/x) + pi/2
     assert asec(x).rewrite(asin) == -asin(1/x) + pi/2
@@ -1755,6 +1780,14 @@ def test_asec_is_real():
     assert asec(t).is_real is False
 
 
+def test_asec_leading_term():
+    assert asec(x).as_leading_term(x, cdir=1) == -I*log(x) + I*log(2)
+    assert asec(x).as_leading_term(x, cdir=-1) == I*log(x) + 2*pi - I*log(2)
+    assert asec(x + 1).as_leading_term(x) == sqrt(2)*sqrt(x)
+    assert asec(x - 1).as_leading_term(x) == pi
+    assert asec(1/x).as_leading_term(x) == pi/2
+
+
 def test_acsc():
     assert acsc(nan) is nan
     assert acsc(1) == pi/2
@@ -1781,7 +1814,6 @@ def test_acsc():
     assert acsc(sqrt(2) - sqrt(6)) == pi*Rational(-5, 12)
 
     assert acsc(x).diff(x) == -1/(x**2*sqrt(1 - 1/x**2))
-    assert acsc(x).as_leading_term(x) == I*log(x)
 
     assert acsc(x).rewrite(log) == -I*log(sqrt(1 - 1/x**2) + I/x)
     assert acsc(x).rewrite(asin) == asin(1/x)
@@ -1808,6 +1840,13 @@ def test_csc_rewrite():
            -1/cos(-pi/2 - 1 + cos(I*besselj(I, I)) +
                   I*cos(-pi/2 + I*besselj(I, I), evaluate=False), evaluate=False)
 
+def test_acsc_leading_term():
+    assert acsc(x).as_leading_term(x, cdir=1) == I*log(x) + pi/2 - I*log(2)
+    assert acsc(x).as_leading_term(x, cdir=-1) == -I*log(x) - 3*pi/2 + I*log(2)
+    assert acsc(x + 1).as_leading_term(x) == pi/2
+    assert acsc(x - 1).as_leading_term(x) == -pi/2
+    assert acsc(1/x).as_leading_term(x) == x
+
 
 def test_inverses_nseries():
     assert asin(x + 2)._eval_nseries(x, 4, None, I) == -asin(2) + pi + sqrt(3)*I*x/3 - sqrt(3)*I*x**2/9 + \
diff --git a/sympy/functions/special/tests/test_bessel.py b/sympy/functions/special/tests/test_bessel.py
index d9ef86d5d7..b0d3089334 100644
--- a/sympy/functions/special/tests/test_bessel.py
+++ b/sympy/functions/special/tests/test_bessel.py
@@ -54,9 +54,9 @@ def test_besselj_leading_term():
 
 
 def test_bessely_leading_term():
-    assert bessely(0, x).as_leading_term(x) == (2*log(x) - 2*log(2))/pi
-    assert bessely(1, sin(x)).as_leading_term(x) == (x*log(x) - x*log(2))/pi
-    assert bessely(1, 2*sqrt(x)).as_leading_term(x) == sqrt(x)*log(x)/pi
+    assert bessely(0, x).as_leading_term(x) == (2*log(x) - 2*log(2) + 2*S.EulerGamma)/pi
+    assert bessely(1, sin(x)).as_leading_term(x) == -2/(pi*x)
+    assert bessely(1, 2*sqrt(x)).as_leading_term(x) == -1/(pi*sqrt(x))
 
 
 def test_besselj_series():
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 3b78c70a25..299109e52d 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -16,7 +16,7 @@
 from sympy.functions.elementary.miscellaneous import (Max, Min, sqrt)
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.trigonometric import (acos, asin, atan, cos, sin, sinc, tan, sec)
-from sympy.functions.special.delta_functions import DiracDelta
+from sympy.functions.special.delta_functions import DiracDelta, Heaviside
 from sympy.functions.special.error_functions import (Ci, Ei, Si, erf, erfc, erfi, fresnelc, li)
 from sympy.functions.special.gamma_functions import (gamma, polygamma)
 from sympy.functions.special.hyper import (hyper, meijerg)
@@ -1724,7 +1724,7 @@ def test_issue_17473():
 def test_issue_17671():
     assert integrate(log(log(x)) / x**2, [x, 1, oo]) == -EulerGamma
     assert integrate(log(log(x)) / x**3, [x, 1, oo]) == -log(2)/2 - EulerGamma/2
-    assert integrate(log(log(x)) / x**10, [x, 1, oo]) == -2*log(3)/9 - EulerGamma/9
+    assert integrate(log(log(x)) / x**10, [x, 1, oo]) == -log(9)/9 - EulerGamma/9
 
 
 def test_issue_2975():
@@ -1892,6 +1892,12 @@ def test_issue_21024():
     assert F_true == integrate(f, x)
 
 
+def test_issue_21721():
+    a = Symbol('a')
+    assert integrate(1/(pi*(1+(x-a)**2)),(x,-oo,oo)).expand() == \
+    -Heaviside(im(a) - 1, 0) + Heaviside(im(a) + 1, 0)
+
+
 def test_issue_21831():
     theta = symbols('theta')
     assert integrate(cos(3*theta)/(5-4*cos(theta)), (theta, 0, 2*pi)) == pi/12
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 0f49719070..656f89149c 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -20,7 +20,7 @@
 from sympy.functions.special.gamma_functions import (digamma, gamma, uppergamma)
 from sympy.integrals.integrals import (Integral, integrate)
 from sympy.series.limits import (Limit, limit)
-from sympy.simplify.simplify import simplify
+from sympy.simplify.simplify import (logcombine, simplify)
 
 from sympy.calculus.accumulationbounds import AccumBounds
 from sympy.core.mul import Mul
@@ -606,6 +606,10 @@ def test_issue_6364():
     assert limit(e, z, 0) == 1/(cos(a)**2 - S.Half)
 
 
+def test_issue_6682():
+    assert limit(exp(2*Ei(-x))/x**2, x, 0) == exp(2*EulerGamma)
+
+
 def test_issue_4099():
     a = Symbol('a')
     assert limit(a/x, x, 0) == oo*sign(a)
@@ -1078,6 +1082,31 @@ def test_issue_20578():
     assert limit(expr,x,0,'+-') == 0
 
 
+def test_issue_21227():
+    f = log(x)
+
+    assert f.nseries(x, logx=y) == y
+    assert f.nseries(x, logx=-x) == -x
+
+    f = log(-log(x))
+
+    assert f.nseries(x, logx=y) == log(-y)
+    assert f.nseries(x, logx=-x) == log(x)
+
+    f = log(log(x))
+
+    assert f.nseries(x, logx=y) == log(y)
+    assert f.nseries(x, logx=-x) == log(-x)
+    assert f.nseries(x, logx=x) == log(x)
+
+    f = log(log(log(1/x)))
+
+    assert f.nseries(x, logx=y) == log(log(-y))
+    assert f.nseries(x, logx=-y) == log(log(y))
+    assert f.nseries(x, logx=x) == log(log(-x))
+    assert f.nseries(x, logx=-x) == log(log(x))
+
+
 def test_issue_21415():
     exp = (x-1)*cos(1/(x-1))
     assert exp.limit(x,1) == 0
@@ -1123,6 +1152,32 @@ def test_issue_22181():
     assert limit((-1)**x * 2**(-x), x, oo) == 0
 
 
+def test_issue_22220():
+    e1 = sqrt(30)*atan(sqrt(30)*tan(x/2)/6)/30
+    e2 = sqrt(30)*I*(-log(sqrt(2)*tan(x/2) - 2*sqrt(15)*I/5) +
+                     +log(sqrt(2)*tan(x/2) + 2*sqrt(15)*I/5))/60
+
+    assert limit(e1, x, -pi) == -sqrt(30)*pi/60
+    assert limit(e2, x, -pi) == -sqrt(30)*pi/30
+
+    assert limit(e1, x, -pi, '-') == sqrt(30)*pi/60
+    assert limit(e2, x, -pi, '-') == 0
+
+    # test https://github.com/sympy/sympy/issues/22220#issuecomment-972727694
+    expr = log(x - I) - log(-x - I)
+    expr2 = logcombine(expr, force=True)
+    assert limit(expr, x, oo) == limit(expr2, x, oo) == I*pi
+
+    # test https://github.com/sympy/sympy/issues/22220#issuecomment-1077618340
+    expr = expr = (-log(tan(x/2) - I) +log(tan(x/2) + I))
+    assert limit(expr, x, pi, '+') == 2*I*pi
+    assert limit(expr, x, pi, '-') == 0
+
+
 def test_issue_23231():
     f = (2**x - 2**(-x))/(2**x + 2**(-x))
     assert limit(f, x, -oo) == -1
+
+
+def test_issue_23596():
+    assert integrate(((1 + x)/x**2)*exp(-1/x), (x, 0, oo)) == oo
diff --git a/sympy/series/tests/test_order.py b/sympy/series/tests/test_order.py
index c3a695f813..d6adaad41e 100644
--- a/sympy/series/tests/test_order.py
+++ b/sympy/series/tests/test_order.py
@@ -157,7 +157,11 @@ def test_add_1():
     assert Order(3*x - 2*x**2) == Order(x)
     assert Order(1 + x) == Order(1, x)
     assert Order(1 + 1/x) == Order(1/x)
-    assert Order(log(x) + 1/log(x)) == Order(log(x))
+    # TODO : A better output for Order(log(x) + 1/log(x))
+    # could be Order(log(x)). Currently Order for expressions
+    # where all arguments would involve a log term would fall
+    # in this category and outputs for these should be improved.
+    assert Order(log(x) + 1/log(x)) == Order((log(x)**2 + 1)/log(x))
     assert Order(exp(1/x) + x) == Order(exp(1/x))
     assert Order(exp(1/x) + 1/x**20) == Order(exp(1/x))
 
@@ -165,7 +169,7 @@ def test_add_1():
 def test_ln_args():
     assert O(log(x)) + O(log(2*x)) == O(log(x))
     assert O(log(x)) + O(log(x**3)) == O(log(x))
-    assert O(log(x*y)) + O(log(x) + log(y)) == O(log(x*y))
+    assert O(log(x*y)) + O(log(x) + log(y)) == O(log(x) + log(y), x, y)
 
 
 def test_multivar_0():