diff --git a/diofant/core/mul.py b/diofant/core/mul.py
index aead0dff87..18901e418d 100644
--- a/diofant/core/mul.py
+++ b/diofant/core/mul.py
@@ -1382,7 +1382,9 @@ def ndiv(a, b):
     def _eval_nseries(self, x, n, logx):
         from ..simplify import powsimp
         terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
-        return powsimp(self.func(*terms).expand(), combine='exp', deep=True)
+        return powsimp(self.func(*terms).expand(power_base=False, power_exp=False,
+                                                mul=False).expand(deep=False),
+                       combine='exp', deep=True)
 
     def _eval_as_leading_term(self, x):
         return self.func(*[t.as_leading_term(x) for t in self.args])
diff --git a/diofant/tests/concrete/test_gosper.py b/diofant/tests/concrete/test_gosper.py
index fffd7d722a..b06e1ac856 100644
--- a/diofant/tests/concrete/test_gosper.py
+++ b/diofant/tests/concrete/test_gosper.py
@@ -50,7 +50,7 @@ def test_gosper_sum():
     # issue sympy/sympy#6033:
     assert gosper_sum(
         n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b)),
-        (n, 0, m)).rewrite(factorial) == \
+        (n, 0, m)).rewrite(factorial).powsimp() == \
         -a*b*(a**m*b**m*factorial(a) *
               factorial(b) - factorial(a + m)*factorial(b + m))/(factorial(a) *
                                                                  factorial(b)*factorial(a + m)*factorial(b + m))
diff --git a/diofant/tests/series/test_limits.py b/diofant/tests/series/test_limits.py
index 26fd6b0feb..8d918c0243 100644
--- a/diofant/tests/series/test_limits.py
+++ b/diofant/tests/series/test_limits.py
@@ -509,7 +509,7 @@ def test_sympyissue_11526():
     df = diff(1/(a*log((x - b)/(x - c))), x)
     res = -1/(-a*c + a*b)
     assert limit(df, x, oo) == res
-    assert limit(simplify(df), x, oo) == res
+    assert (limit(simplify(df), x, oo) - res).simplify() == 0
 
     e = log((1/x - b)/(1/x - c))
     assert e.as_leading_term(x) == x*(c - b)
diff --git a/diofant/tests/series/test_residues.py b/diofant/tests/series/test_residues.py
index 2853bf9814..c8595c604b 100644
--- a/diofant/tests/series/test_residues.py
+++ b/diofant/tests/series/test_residues.py
@@ -1,7 +1,7 @@
 import pytest
 
-from diofant import (Function, I, Rational, Symbol, cot, exp, factorial, log,
-                     pi, residue, root, sin, sqrt, tanh)
+from diofant import (Function, I, Mul, Rational, Symbol, cot, exp, factorial,
+                     log, pi, residue, root, sin, sqrt, tan, tanh)
 from diofant.abc import a, s, x, z
 
 
@@ -85,3 +85,11 @@ def test_sympyissue_21177():
 
     assert residue(e1, x, pt) == ans
     assert residue(e2, x, pt) == ans
+
+
+def test_sympyissue_21176():
+    e = x**2*cot(pi*x)/(x**4 + 1)
+    pt = -sqrt(2)/2 - sqrt(2)*I/2
+    assert residue(e, x, pt) == sqrt(2)*I/Mul(2, -2 + 2*I, tan(sqrt(2)*pi/2 +
+                                                               sqrt(2)*I*pi/2),
+                                              evaluate=False)
diff --git a/diofant/tests/series/test_series.py b/diofant/tests/series/test_series.py
index b0aa43b5ba..80c9108e2a 100644
--- a/diofant/tests/series/test_series.py
+++ b/diofant/tests/series/test_series.py
@@ -261,12 +261,12 @@ def test_sympyissue_21245():
     e = 1/(1 - x - x**2)
     assert (e.series(x, 1/fi, 2) ==
             -sqrt(5)/(Mul(5, x - 1/(Rational(1, 2) + sqrt(5)/2),
-                      evaluate=False)) + sqrt(5)/(5 + 5*sqrt(5)) +
-            1/(5 + 5*sqrt(5)) +
-            (x - 1/(Rational(1, 2) + sqrt(5)/2))*(-6*sqrt(5)/(50*sqrt(5) + 150) -
-                                                  10/(50*sqrt(5) + 150)) +
-            O((x - sqrt(5)/2 + Rational(1, 2))**2,
-              (x, -Rational(1, 2) + sqrt(5)/2)))
+                          evaluate=False)) + 1/(sqrt(5) + 5) +
+            sqrt(5)/(Mul(5, sqrt(5) + 5, evaluate=False)) +
+            (x - 1/(Rational(1, 2) + sqrt(5)/2)) *
+            (-6*sqrt(5)/(Mul(5, 10*sqrt(5) + 30, evaluate=False)) -
+             2/(10*sqrt(5) + 30)) + O((x - sqrt(5)/2 + Rational(1, 2))**2,
+                                      (x, -Rational(1, 2) + sqrt(5)/2)))
 
 
 def test_diofantissue_1139():
diff --git a/diofant/tests/test_wester.py b/diofant/tests/test_wester.py
index 835c18c3d7..6a6ba7052e 100644
--- a/diofant/tests/test_wester.py
+++ b/diofant/tests/test_wester.py
@@ -2017,7 +2017,8 @@ def test_X6():
     # Taylor series of nonscalar objects (noncommutative multiplication)
     # expected result => (B A - A B) t^2/2 + O(t^3)   [Stanly Steinberg]
     a, b = symbols('a b', commutative=False, scalar=False)
-    assert (series(exp((a + b)*x) - exp(a*x) * exp(b*x), x, x0=0, n=3) ==
+    assert (series(exp((a + b)*x) - exp(a*x) * exp(b*x), x,
+                   x0=0, n=3).expand().collect(x) ==
             x**2*(-a*b/2 + b*a/2) + O(x**3))
 
 
diff --git a/docs/release/notes-0.13.rst b/docs/release/notes-0.13.rst
index 6af68c2471..a5c3358cb9 100644
--- a/docs/release/notes-0.13.rst
+++ b/docs/release/notes-0.13.rst
@@ -116,3 +116,4 @@ These Sympy issues also were addressed:
 * :sympyissue:`21785`: Limit gives TypeError from as_leading_term
 * :sympyissue:`21812`: LambertW displaying in jupyter lab
 * :sympyissue:`21814`: Printing of unevaluated Mul needs brackets
+* :sympyissue:`21176`: Incorrect residue of x**2*cot(pi*x)/(x**4 + 1)