Adds _eval_is_meromorphic() to uppergamma and lowergamma function

sympy/functions/special/bessel.py

@@ -460,12 +460,10 @@ def eval(cls, nu, z):
                 return S.ComplexInfinity
             elif re(nu).is_zero:
                 return S.NaN
-        if im(z) is S.Infinity or im(z) is S.NegativeInfinity:
+        if z in (S.Infinity, I*S.Infinity, I*S.NegativeInfinity):
             return S.Zero
 
         if nu.is_integer:
-            if z is S.Infinity:
-                return S.Zero
             if nu.could_extract_minus_sign():
                 return besselk(-nu, z)
 

sympy/functions/special/gamma_functions.py

@@ -364,13 +364,30 @@ def _eval_conjugate(self):
         if x not in (S.Zero, S.NegativeInfinity):
             return self.func(self.args[0].conjugate(), x.conjugate())
 
+    def _eval_is_meromorphic(self, x, a):
+        # By https://en.wikipedia.org/wiki/Incomplete_gamma_function#Holomorphic_extension,
+        #    lowergamma(s, z) = z**s*gamma(s)*gammastar(s, z),
+        # where gammastar(s, z) is holomorphic for all s and z.
+        # Hence the singularities of lowergamma are z = 0  (branch
+        # point) and nonpositive integer values of s (poles of gamma(s)).
+        s, z = self.args
+        args_merom = fuzzy_and([z._eval_is_meromorphic(x, a),
+            s._eval_is_meromorphic(x, a)])
+        if not args_merom:
+            return args_merom
+        z0 = z.subs(x, a)
+        if s.is_integer:
+            return fuzzy_and([s.is_positive, z0.is_finite])
+        s0 = s.subs(x, a)
+        return fuzzy_and([s0.is_finite, z0.is_finite, fuzzy_not(z0.is_zero)])
+
     def _eval_aseries(self, n, args0, x, logx):
         from sympy import O
         s, z = self.args
         if args0[0] is S.Infinity and not z.has(x):
             coeff = z**s*exp(-z)
             sum_expr = sum(z**k/rf(s, k + 1) for k in range(n - 1))
-            o = O(s**(-n - 1))
+            o = O(z**s*s**(-n))
             return coeff*sum_expr + o
         return super()._eval_aseries(n, args0, x, logx)
 
@@ -533,6 +550,9 @@ def _eval_conjugate(self):
         if not z in (S.Zero, S.NegativeInfinity):
             return self.func(self.args[0].conjugate(), z.conjugate())
 
+    def _eval_is_meromorphic(self, x, a):
+        return lowergamma._eval_is_meromorphic(self, x, a)
+
     def _eval_rewrite_as_lowergamma(self, s, x, **kwargs):
         return gamma(s) - lowergamma(s, x)
 

sympy/functions/special/tests/test_gamma_functions.py

@@ -140,8 +140,27 @@ def test_lowergamma():
     assert conjugate(lowergamma(x, 0)) == 0
     assert unchanged(conjugate, lowergamma(x, -oo))
 
-    assert lowergamma(x, 1).series(x, oo, 3) == \
-        (1 + 1/(x + 1))*exp(-1)/x + O(x**(-3), (x, oo))
+    assert lowergamma(0, x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(S(1)/3, x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(1, x, evaluate=False)._eval_is_meromorphic(x, 0) == True
+    assert lowergamma(x, x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(x + 1, x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(1/x, x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(0, x + 1)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(S(1)/3, x + 1)._eval_is_meromorphic(x, 0) == True
+    assert lowergamma(1, x + 1, evaluate=False)._eval_is_meromorphic(x, 0) == True
+    assert lowergamma(x, x + 1)._eval_is_meromorphic(x, 0) == True
+    assert lowergamma(x + 1, x + 1)._eval_is_meromorphic(x, 0) == True
+    assert lowergamma(1/x, x + 1)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(0, 1/x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(S(1)/3, 1/x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(1, 1/x, evaluate=False)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(x, 1/x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(x + 1, 1/x)._eval_is_meromorphic(x, 0) == False
+    assert lowergamma(1/x, 1/x)._eval_is_meromorphic(x, 0) == False
+
+    assert lowergamma(x, 2).series(x, oo, 3) == \
+        2**x*(1 + 2/(x + 1))*exp(-2)/x + O(exp(x*log(2))/x**3, (x, oo))
 
     assert lowergamma(
         x, y).rewrite(expint) == -y**x*expint(-x + 1, y) + gamma(x)