Repository: sympy/sympy

Issue:
`sin(Add(1, -1, evaluate=False)).is_zero` returns `False`
Should return `None` or ideally `True`.

Note that `sin(Add(pi/4, pi/4, evaluate=False)).is_zero` gives the correct answer (`True`), but also that `sin(Add(pi/4, pi/4, evaluate=False))` actually evaluates to `0`. Hence, it seems to be that for `Add(1, -1, evaluate=False)`, `sin` does not evaluate.

Compare this to `sqrt(Add(1, -1, evaluate=False))` which doesn't evaluate, but returns `True` for `is_zero`.

So there is both an inconsistency in how `sin` and `sqrt` evaluates, but more importantly, the `is_zero` result from `sin` is incorrect.

As stated in #17510 this can be easily solved by adding an `is_zero` method for `sin`. However, this will only hide the underlying reason of the problem which is likely to pop up for some other function, so a fix should not include adding `_eval_is_zero` for `sin` (unless there are very good arguments for that the problem is related to `sin` and only `sin`).

Also, note that `Add(1, -1, evaluate=False).is_zero` is `True`, so there is something with the combination of (at least) `sin` and this expression.

(Furthermore, there is an XFAIL-test for this in #17510.)

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v147). 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:
* core
  * `is_zero` and `is_rational` more often returns a correct answer for unevaluated expressions. ([#17510](https://github.com/sympy/sympy/pull/17510) by [@oscargus](https://github.com/oscargus))

* functions
  * `hyper` and `appellf1` supports `evaluate` keyword. ([#17510](https://github.com/sympy/sympy/pull/17510) by [@oscargus](https://github.com/oscargus))

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

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    …appelf1 and hyper

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234". See
    https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
    write a comment on that issue linking back to this pull request once it is
    open. -->


    #### Brief description of what is fixed or changed
    Now, expressions as `sin(0, evaluate=False).is_zero` provides a correct answer (holds for all function). Added something similar for `is_rational`.

    As part of it, I found an issue `sin(Add(1, .1, evaluate=False)).is_zero` returns `False`, but not the cause of this PR. Added a test for it as well. (Can be solved by adding a specific `is_zero` function to `sin`, but seems like a problem that will just pop up in another function etc...)

    Also, added support for `evaluate`-keywords to `hyper` and `appellf1`. `hyper` never evaluates, but I think it makes sense to support it anyway, if nothing else for the `is_zero`, ìs_rational`-support.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. See
    https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
    on how to write release notes. The bot will check your release notes
    automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    * core
        * `is_zero` and `is_rational` more often returns a correct answer for unevaluated expressions.

    * functions
       * `hyper` and `appellf1` supports `evaluate` keyword.       
    <!-- END RELEASE NOTES -->


</details><p>

oscargus: Turns out that evaluating huge factorials isn't really feasible... So I'll have to rethink this a bit...
oscarbenjamin: > Turns out that evaluating huge factorials isn't really feasible... So I'll have to rethink this a bit...

This is what I meant when I said that the user probably has a good reason for using evaluate=False :).

As a general principle I think that the assumptions system should lie underneath the evaluation system. Evaluation should use the assumptions system to make decisions but the assumptions system should not use evaluation to answer queries.
smichr: > the assumptions system should not use evaluation to answer queries

If this is the case then how should a query about whether some complicated expression is positive be handled when `expr.is_positive is None`? Maybe that is where SymPy should recommend that users use `expr > 0` to get evaluation, if possible (since it might not be possible if `expr` cannot be evaluated with precision).
oscargus: I can see the point of not using the evaluation system in `_eval_is_xx` methods. However, it is worth noting that it currently is doing that in some places.

> If this is the case then how should a query about whether some complicated expression is positive be handled when expr.is_positive is None? Maybe that is where SymPy should recommend that users use expr > 0 to get evaluation, if possible (since it might not be possible if expr cannot be evaluated with precision).

I think that in the code base one can see all types of ways to solve this and I think that one should be the preferred.

Btw, this PR has now changed focus to implementing better `eval` and `_eval_is_zero` for functions.

> Given that in a parallel thread you are talking about mechanically changing the code to things that are faster I should point out the main difference here which is that is_zero can take an arbitrary amount of time.

As you see, in this PR I (now) first check `is S.Zero` (which is faster than `== 0` and equivalent as long as the functions returns SymPy objects which they should...) and then later with `is_zero`. Except that I haven't rewritten cases that only checks with `is_zero`. If one can detect a zero valued expression early on that should be beneficial from an execution time perspective compared to propagating symbolic expressions.
oscargus: Btw, I'll add tests and revert the `quantum` fix. Although it should be relevant for consistency to be able to provide the `evaluate` keyword to `Dagger`, I cannot really figure out the problem, plus that the printing assumes that it is always evaluated. So not really worth the hassle right now.
oscarbenjamin: > I can see the point of not using the evaluation system in `_eval_is_xx` methods. However, it is worth noting that it currently is doing that in some places.

Can you show an example?

There are cases where a subexpression is created for example `if (self.args[0]/pi).is_integer:` Is that what you mean?
oscargus: Search for `self.func(*self.args)` to see all cases. Some are potentially problematic, although none enforces `evaluate=True`.

Seems like `log` is the worst case when it comes to `_is_eval_xx`, e.g.
```
    def _eval_is_rational(self):
        s = self.func(*self.args)
        if s.func == self.func:
            if (self.args[0] - 1).is_zero:
                return True
            if s.args[0].is_rational and fuzzy_not((self.args[0] - 1).is_zero):
                return False
        else:
            return s.is_rational
```

But the maybe most worrying is `Basic.copy`:

```
a = Add(1, 1, evaluate=True)
a.copy() # Results in 2
```

Now, I'm not sure what the purpose of copy is, but I can imagine that if someone uses it, there may be unexpected consequences.
oscarbenjamin: > Seems like `log` is the worst case when it comes to `_is_eval_xx`, e.g.
> 
> ```
>     def _eval_is_rational(self):
>         s = self.func(*self.args)
> ```

This should be removed I think.

> Now, I'm not sure what the purpose of copy is, but I can imagine that if someone uses it, there may be unexpected consequences.

There shouldn't be any need for Basic.copy since Basic is immutable. Either way though it is difficult at that level to use an evaluate kwarg because many Basic subclasses don't support it #17280 #17372.
codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/17510?src=pr&el=h1) Report
> Merging [#17510](https://codecov.io/gh/sympy/sympy/pull/17510?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/f88589f4512870794834f57edb39ec1afa730b68?src=pr&el=desc) will **increase** coverage by `7.754%`.
> The diff coverage is `45.49%`.

```diff
@@              Coverage Diff              @@
##            master    #17510       +/-   ##
=============================================
+ Coverage   66.893%   74.648%   +7.754%     
=============================================
  Files          633       633               
  Lines       164469    164667      +198     
  Branches     38616     38697       +81     
=============================================
+ Hits        110019    122921    +12902     
+ Misses       49027     36286    -12741     
- Partials      5423      5460       +37
```

oscarbenjamin: This generally looks good although I left comments.

There are lots of places where we now have something like:
```python
def eval(cls, arg):
    if arg is S.Zero:
        # Some stuff
    # Some other stuff
    if arg.is_zero:
        # Same stuff as before
```
I think it would be better to combine those. I presume that you do this because it will possibly be quicker but I don't think it will be because in most cases where another condition applies the `is_zero` check will be cheap and `S.Zero.is_zero` is also potentially cheaper than `S.Zero is S.Zero`.
oscarbenjamin: @oscargus are you still working on this? If not I think this should be merged.

Patch:
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index bfee2c2e89..faae1b8a93 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -319,6 +319,9 @@ def eval(cls, arg):
         elif isinstance(arg, MatrixBase):
             return arg.exp()
 
+        if arg.is_zero:
+            return S.One
+
     @property
     def base(self):
         """
@@ -642,6 +645,9 @@ def eval(cls, arg, base=None):
             elif arg is S.Exp1:
                 return S.One
 
+        if arg.is_zero:
+            return S.ComplexInfinity
+
         # don't autoexpand Pow or Mul (see the issue 3351):
         if not arg.is_Add:
             coeff = arg.as_coefficient(I)
@@ -962,6 +968,8 @@ def eval(cls, x, k=None):
                 return S.Exp1
             if x is S.Infinity:
                 return S.Infinity
+            if x.is_zero:
+                return S.Zero
 
         if fuzzy_not(k.is_zero):
             if x.is_zero:
@@ -1038,3 +1046,12 @@ def _eval_nseries(self, x, n, logx):
 
             return s + Order(x**n, x)
         return super(LambertW, self)._eval_nseries(x, n, logx)
+
+    def _eval_is_zero(self):
+        x = self.args[0]
+        if len(self.args) == 1:
+            k = S.Zero
+        else:
+            k = self.args[1]
+        if x.is_zero and k.is_zero:
+            return True
diff --git a/sympy/functions/elementary/hyperbolic.py b/sympy/functions/elementary/hyperbolic.py
index 526d83c7fe..fef5464c91 100644
--- a/sympy/functions/elementary/hyperbolic.py
+++ b/sympy/functions/elementary/hyperbolic.py
@@ -127,6 +127,9 @@ def eval(cls, arg):
                 if m:
                     return sinh(m)*cosh(x) + cosh(m)*sinh(x)
 
+            if arg.is_zero:
+                return S.Zero
+
             if arg.func == asinh:
                 return arg.args[0]
 
@@ -243,6 +246,11 @@ def _eval_is_finite(self):
         if arg.is_imaginary:
             return True
 
+    def _eval_is_zero(self):
+        arg = self.args[0]
+        if arg.is_zero:
+            return True
+
 
 class cosh(HyperbolicFunction):
     r"""
@@ -295,6 +303,9 @@ def eval(cls, arg):
                 if m:
                     return cosh(m)*cosh(x) + sinh(m)*sinh(x)
 
+            if arg.is_zero:
+                return S.One
+
             if arg.func == asinh:
                 return sqrt(1 + arg.args[0]**2)
 
@@ -460,6 +471,9 @@ def eval(cls, arg):
                     else: # tanhm == 0
                         return tanh(x)
 
+            if arg.is_zero:
+                return S.Zero
+
             if arg.func == asinh:
                 x = arg.args[0]
                 return x/sqrt(1 + x**2)
@@ -551,6 +565,11 @@ def _eval_is_finite(self):
         if arg.is_extended_real:
             return True
 
+    def _eval_is_zero(self):
+        arg = self.args[0]
+        if arg.is_zero:
+            return True
+
 
 class coth(HyperbolicFunction):
     r"""
@@ -610,6 +629,9 @@ def eval(cls, arg):
                     else: # cothm == 0
                         return tanh(x)
 
+            if arg.is_zero:
+                return S.ComplexInfinity
+
             if arg.func == asinh:
                 x = arg.args[0]
                 return sqrt(1 + x**2)/x
@@ -919,6 +941,9 @@ def eval(cls, arg):
             if arg is S.ComplexInfinity:
                 return S.ComplexInfinity
 
+            if arg.is_zero:
+                return S.Zero
+
             i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
@@ -975,6 +1000,11 @@ def inverse(self, argindex=1):
         """
         return sinh
 
+    def _eval_is_zero(self):
+        arg = self.args[0]
+        if arg.is_zero:
+            return True
+
 
 class acosh(InverseHyperbolicFunction):
     """
@@ -1048,6 +1078,9 @@ def eval(cls, arg):
         if arg == -S.ImaginaryUnit*S.Infinity:
             return S.Infinity - S.ImaginaryUnit*S.Pi/2
 
+        if arg.is_zero:
+            return S.Pi*S.ImaginaryUnit*S.Half
+
         if isinstance(arg, cosh) and arg.args[0].is_number:
             z = arg.args[0]
             if z.is_real:
@@ -1158,6 +1191,9 @@ def eval(cls, arg):
                 if _coeff_isneg(arg):
                     return -cls(-arg)
 
+        if arg.is_zero:
+            return S.Zero
+
         if isinstance(arg, tanh) and arg.args[0].is_number:
             z = arg.args[0]
             if z.is_real:
@@ -1193,6 +1229,12 @@ def _eval_as_leading_term(self, x):
     def _eval_rewrite_as_log(self, x, **kwargs):
         return (log(1 + x) - log(1 - x)) / 2
 
+    def _eval_is_zero(self):
+        arg = self.args[0]
+        if arg.is_zero:
+            return True
+
+
     def inverse(self, argindex=1):
         """
         Returns the inverse of this function.
@@ -1245,6 +1287,9 @@ def eval(cls, arg):
                 if _coeff_isneg(arg):
                     return -cls(-arg)
 
+        if arg.is_zero:
+            return S.Pi*S.ImaginaryUnit*S.Half
+
     @staticmethod
     @cacheit
     def taylor_term(n, x, *previous_terms):
@@ -1375,6 +1420,9 @@ def eval(cls, arg):
             from sympy.calculus.util import AccumBounds
             return S.ImaginaryUnit*AccumBounds(-S.Pi/2, S.Pi/2)
 
+        if arg.is_zero:
+            return S.Infinity
+
     @staticmethod
     @cacheit
     def expansion_term(n, x, *previous_terms):
@@ -1484,6 +1532,9 @@ def eval(cls, arg):
         if arg is S.ComplexInfinity:
             return S.Zero
 
+        if arg.is_zero:
+            return S.ComplexInfinity
+
         if _coeff_isneg(arg):
             return -cls(-arg)
 
diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py
index 179e705557..f1d75e8593 100644
--- a/sympy/functions/elementary/trigonometric.py
+++ b/sympy/functions/elementary/trigonometric.py
@@ -196,6 +196,8 @@ def _pi_coeff(arg, cycles=1):
                 else:
                     return c2*x
             return cx
+    elif arg.is_zero:
+        return S.Zero
 
 
 class sin(TrigonometricFunction):
@@ -339,6 +341,9 @@ def eval(cls, arg):
             if m:
                 return sin(m)*cos(x) + cos(m)*sin(x)
 
+        if arg.is_zero:
+            return S.Zero
+
         if isinstance(arg, asin):
             return arg.args[0]
 
@@ -476,6 +481,11 @@ def _eval_is_finite(self):
         if arg.is_extended_real:
             return True
 
+    def _eval_is_zero(self):
+        arg = self.args[0]
+        if arg.is_zero:
+            return True
+
 
 class cos(TrigonometricFunction):
     """
@@ -643,6 +653,9 @@ def eval(cls, arg):
             if m:
                 return cos(m)*cos(x) - sin(m)*sin(x)
 
+        if arg.is_zero:
+            return S.One
+
         if isinstance(arg, acos):
             return arg.args[0]
 
@@ -1063,6 +1076,9 @@ def eval(cls, arg):
                 else: # tanm == 0
                     return tan(x)
 
+        if arg.is_zero:
+            return S.Zero
+
         if isinstance(arg, atan):
             return arg.args[0]
 
@@ -1222,6 +1238,11 @@ def _eval_is_finite(self):
         if arg.is_imaginary:
             return True
 
+    def _eval_is_zero(self):
+        arg = self.args[0]
+        if arg.is_zero:
+            return True
+
 
 class cot(TrigonometricFunction):
     """
@@ -1356,6 +1377,9 @@ def eval(cls, arg):
                 else: # cotm == 0
                     return -tan(x)
 
+        if arg.is_zero:
+            return S.ComplexInfinity
+
         if isinstance(arg, acot):
             return arg.args[0]
 
@@ -2059,6 +2083,9 @@ def eval(cls, arg):
         if i_coeff is not None:
             return S.ImaginaryUnit * asinh(i_coeff)
 
+        if arg.is_zero:
+            return S.Zero
+
         if isinstance(arg, sin):
             ang = arg.args[0]
             if ang.is_comparable:
@@ -2403,6 +2430,9 @@ def eval(cls, arg):
         if i_coeff is not None:
             return S.ImaginaryUnit * atanh(i_coeff)
 
+        if arg.is_zero:
+            return S.Zero
+
         if isinstance(arg, tan):
             ang = arg.args[0]
             if ang.is_comparable:
@@ -2572,6 +2602,9 @@ def eval(cls, arg):
         if i_coeff is not None:
             return -S.ImaginaryUnit * acoth(i_coeff)
 
+        if arg.is_zero:
+            return S.Pi*S.Half
+
         if isinstance(arg, cot):
             ang = arg.args[0]
             if ang.is_comparable:
diff --git a/sympy/functions/special/bessel.py b/sympy/functions/special/bessel.py
index 3eaa0a2e41..08b3d2bbbd 100644
--- a/sympy/functions/special/bessel.py
+++ b/sympy/functions/special/bessel.py
@@ -1117,6 +1117,8 @@ def eval(cls, arg):
                 return S.Zero
             elif arg.is_zero:
                 return S.One / (3**Rational(2, 3) * gamma(Rational(2, 3)))
+        if arg.is_zero:
+            return S.One / (3**Rational(2, 3) * gamma(Rational(2, 3)))
 
     def fdiff(self, argindex=1):
         if argindex == 1:
@@ -1286,6 +1288,9 @@ def eval(cls, arg):
             elif arg.is_zero:
                 return S.One / (3**Rational(1, 6) * gamma(Rational(2, 3)))
 
+        if arg.is_zero:
+            return S.One / (3**Rational(1, 6) * gamma(Rational(2, 3)))
+
     def fdiff(self, argindex=1):
         if argindex == 1:
             return airybiprime(self.args[0])
@@ -1440,8 +1445,9 @@ def eval(cls, arg):
                 return S.NaN
             elif arg is S.Infinity:
                 return S.Zero
-            elif arg.is_zero:
-                return -S.One / (3**Rational(1, 3) * gamma(Rational(1, 3)))
+
+        if arg.is_zero:
+            return S.NegativeOne / (3**Rational(1, 3) * gamma(Rational(1, 3)))
 
     def fdiff(self, argindex=1):
         if argindex == 1:
@@ -1599,6 +1605,10 @@ def eval(cls, arg):
             elif arg.is_zero:
                 return 3**Rational(1, 6) / gamma(Rational(1, 3))
 
+        if arg.is_zero:
+            return 3**Rational(1, 6) / gamma(Rational(1, 3))
+
+
     def fdiff(self, argindex=1):
         if argindex == 1:
             return self.args[0]*airybi(self.args[0])
@@ -1710,20 +1720,28 @@ class marcumq(Function):
     @classmethod
     def eval(cls, m, a, b):
         from sympy import exp, uppergamma
-        if a == 0:
-            if m == 0 and b == 0:
+        if a is S.Zero:
+            if m is S.Zero and b is S.Zero:
                 return S.Zero
-            return uppergamma(m, b**2 / 2) / gamma(m)
+            return uppergamma(m, b**2 * S.Half) / gamma(m)
 
-        if m == 0 and b == 0:
-            return 1 - 1 / exp(a**2 / 2)
+        if m is S.Zero and b is S.Zero:
+            return 1 - 1 / exp(a**2 * S.Half)
 
         if a == b:
-            if m == 1:
-                return (1 + exp(-a**2) * besseli(0, a**2)) / 2
+            if m is S.One:
+                return (1 + exp(-a**2) * besseli(0, a**2))*S.Half
             if m == 2:
                 return S.Half + S.Half * exp(-a**2) * besseli(0, a**2) + exp(-a**2) * besseli(1, a**2)
 
+        if a.is_zero:
+            if m.is_zero and b.is_zero:
+                return S.Zero
+            return uppergamma(m, b**2*S.Half) / gamma(m)
+
+        if m.is_zero and b.is_zero:
+            return 1 - 1 / exp(a**2*S.Half)
+
     def fdiff(self, argindex=2):
         from sympy import exp
         m, a, b = self.args
@@ -1753,3 +1771,7 @@ def _eval_rewrite_as_besseli(self, m, a, b, **kwargs):
             if m.is_Integer and m >= 2:
                 s = sum([besseli(i, a**2) for i in range(1, m)])
                 return S.Half + exp(-a**2) * besseli(0, a**2) / 2 + exp(-a**2) * s
+
+    def _eval_is_zero(self):
+        if all(arg.is_zero for arg in self.args):
+            return True
diff --git a/sympy/functions/special/beta_functions.py b/sympy/functions/special/beta_functions.py
index 2d78bd25e5..fda042bc3a 100644
--- a/sympy/functions/special/beta_functions.py
+++ b/sympy/functions/special/beta_functions.py
@@ -1,5 +1,6 @@
 from __future__ import print_function, division
 
+from sympy.core import S
 from sympy.core.function import Function, ArgumentIndexError
 from sympy.functions.special.gamma_functions import gamma, digamma
 
@@ -33,35 +34,35 @@ class beta(Function):
     ========
 
     >>> from sympy import I, pi
-    >>> from sympy.abc import x,y
+    >>> from sympy.abc import x, y
 
     The Beta function obeys the mirror symmetry:
 
     >>> from sympy import beta
     >>> from sympy import conjugate
-    >>> conjugate(beta(x,y))
+    >>> conjugate(beta(x, y))
     beta(conjugate(x), conjugate(y))
 
     Differentiation with respect to both x and y is supported:
 
     >>> from sympy import beta
     >>> from sympy import diff
-    >>> diff(beta(x,y), x)
+    >>> diff(beta(x, y), x)
     (polygamma(0, x) - polygamma(0, x + y))*beta(x, y)
 
     >>> from sympy import beta
     >>> from sympy import diff
-    >>> diff(beta(x,y), y)
+    >>> diff(beta(x, y), y)
     (polygamma(0, y) - polygamma(0, x + y))*beta(x, y)
 
     We can numerically evaluate the gamma function to arbitrary precision
     on the whole complex plane:
 
     >>> from sympy import beta
-    >>> beta(pi,pi).evalf(40)
+    >>> beta(pi, pi).evalf(40)
     0.02671848900111377452242355235388489324562
 
-    >>> beta(1+I,1+I).evalf(20)
+    >>> beta(1 + I, 1 + I).evalf(20)
     -0.2112723729365330143 - 0.7655283165378005676*I
 
     See Also
@@ -98,7 +99,10 @@ def fdiff(self, argindex):
 
     @classmethod
     def eval(cls, x, y):
-        pass
+        if y is S.One:
+            return 1/x
+        if x is S.One:
+            return 1/y
 
     def _eval_expand_func(self, **hints):
         x, y = self.args
diff --git a/sympy/functions/special/delta_functions.py b/sympy/functions/special/delta_functions.py
index b592b9a395..e2388400a8 100644
--- a/sympy/functions/special/delta_functions.py
+++ b/sympy/functions/special/delta_functions.py
@@ -202,12 +202,12 @@ def eval(cls, arg, k=0):
                 Complex part: %s  found in %s .''' % (
                 repr(im(arg)), repr(arg))))
         c, nc = arg.args_cnc()
-        if c and c[0] == -1:
+        if c and c[0] is S.NegativeOne:
             # keep this fast and simple instead of using
             # could_extract_minus_sign
-            if k % 2 == 1:
+            if k.is_odd:
                 return -cls(-arg, k)
-            elif k % 2 == 0:
+            elif k.is_even:
                 return cls(-arg, k) if k else cls(-arg)
 
     @deprecated(useinstead="expand(diracdelta=True, wrt=x)", issue=12859, deprecated_since_version="1.1")
diff --git a/sympy/functions/special/elliptic_integrals.py b/sympy/functions/special/elliptic_integrals.py
index 978a8b4fe7..adee570a6c 100644
--- a/sympy/functions/special/elliptic_integrals.py
+++ b/sympy/functions/special/elliptic_integrals.py
@@ -67,6 +67,9 @@ def eval(cls, m):
                    I*S.NegativeInfinity, S.ComplexInfinity):
             return S.Zero
 
+        if m.is_zero:
+            return pi*S.Half
+
     def fdiff(self, argindex=1):
         m = self.args[0]
         return (elliptic_e(m) - (1 - m)*elliptic_k(m))/(2*m*(1 - m))
@@ -81,11 +84,16 @@ def _eval_nseries(self, x, n, logx):
         return hyperexpand(self.rewrite(hyper)._eval_nseries(x, n=n, logx=logx))
 
     def _eval_rewrite_as_hyper(self, m, **kwargs):
-        return (pi/2)*hyper((S.Half, S.Half), (S.One,), m)
+        return pi*S.Half*hyper((S.Half, S.Half), (S.One,), m)
 
     def _eval_rewrite_as_meijerg(self, m, **kwargs):
         return meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -m)/2
 
+    def _eval_is_zero(self):
+        m = self.args[0]
+        if m.is_infinite:
+            return True
+
     def _sage_(self):
         import sage.all as sage
         return sage.elliptic_kc(self.args[0]._sage_())
@@ -132,12 +140,12 @@ class elliptic_f(Function):
 
     @classmethod
     def eval(cls, z, m):
-        k = 2*z/pi
+        if z.is_zero:
+            return S.Zero
         if m.is_zero:
             return z
-        elif z.is_zero:
-            return S.Zero
-        elif k.is_integer:
+        k = 2*z/pi
+        if k.is_integer:
             return k*elliptic_k(m)
         elif m in (S.Infinity, S.NegativeInfinity):
             return S.Zero
@@ -159,6 +167,13 @@ def _eval_conjugate(self):
         if (m.is_real and (m - 1).is_positive) is False:
             return self.func(z.conjugate(), m.conjugate())
 
+    def _eval_is_zero(self):
+        z, m = self.args
+        if z.is_zero:
+            return True
+        if m.is_extended_real and m.is_infinite:
+            return True
+
 
 class elliptic_e(Function):
     r"""
@@ -319,14 +334,14 @@ class elliptic_pi(Function):
     def eval(cls, n, m, z=None):
         if z is not None:
             n, z, m = n, m, z
-            k = 2*z/pi
             if n.is_zero:
                 return elliptic_f(z, m)
-            elif n == S.One:
+            elif n is S.One:
                 return (elliptic_f(z, m) +
                         (sqrt(1 - m*sin(z)**2)*tan(z) -
                          elliptic_e(z, m))/(1 - m))
-            elif k.is_integer:
+            k = 2*z/pi
+            if k.is_integer:
                 return k*elliptic_pi(n, m)
             elif m.is_zero:
                 return atanh(sqrt(n - 1)*tan(z))/sqrt(n - 1)
@@ -339,10 +354,15 @@ def eval(cls, n, m, z=None):
                 return S.Zero
             elif z.could_extract_minus_sign():
                 return -elliptic_pi(n, -z, m)
+            if n.is_zero:
+                return elliptic_f(z, m)
+            if m.is_extended_real and m.is_infinite or \
+                    n.is_extended_real and n.is_infinite:
+                return S.Zero
         else:
             if n.is_zero:
                 return elliptic_k(m)
-            elif n == S.One:
+            elif n is S.One:
                 return S.ComplexInfinity
             elif m.is_zero:
                 return pi/(2*sqrt(1 - n))
@@ -354,6 +374,11 @@ def eval(cls, n, m, z=None):
                 return S.Zero
             elif m in (S.Infinity, S.NegativeInfinity):
                 return S.Zero
+            if n.is_zero:
+                return elliptic_k(m)
+            if m.is_extended_real and m.is_infinite or \
+                    n.is_extended_real and n.is_infinite:
+                return S.Zero
 
     def _eval_conjugate(self):
         if len(self.args) == 3:
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 2d306fb78f..220de4eb36 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -145,6 +145,9 @@ def eval(cls, arg):
         if isinstance(arg, erfcinv):
             return S.One - arg.args[0]
 
+        if arg.is_zero:
+            return S.Zero
+
         # Only happens with unevaluated erf2inv
         if isinstance(arg, erf2inv) and arg.args[0].is_zero:
             return arg.args[1]
@@ -183,6 +186,10 @@ def _eval_is_finite(self):
         else:
             return self.args[0].is_extended_real
 
+    def _eval_is_zero(self):
+        if self.args[0].is_zero:
+            return True
+
     def _eval_rewrite_as_uppergamma(self, z, **kwargs):
         from sympy import uppergamma
         return sqrt(z**2)/z*(S.One - uppergamma(S.Half, z**2)/sqrt(S.Pi))
@@ -326,6 +333,9 @@ def eval(cls, arg):
         if isinstance(arg, erfcinv):
             return arg.args[0]
 
+        if arg.is_zero:
+            return S.One
+
         # Try to pull out factors of I
         t = arg.extract_multiplicatively(S.ImaginaryUnit)
         if t is S.Infinity or t is S.NegativeInfinity:
@@ -487,6 +497,9 @@ def eval(cls, z):
             elif z is S.Infinity:
                 return S.Infinity
 
+        if z.is_zero:
+            return S.Zero
+
         # Try to pull out factors of -1
         if z.could_extract_minus_sign():
             return -cls(-z)
@@ -523,6 +536,10 @@ def _eval_conjugate(self):
     def _eval_is_extended_real(self):
         return self.args[0].is_extended_real
 
+    def _eval_is_zero(self):
+        if self.args[0].is_zero:
+            return True
+
     def _eval_rewrite_as_tractable(self, z, **kwargs):
         return self.rewrite(erf).rewrite("tractable", deep=True)
 
@@ -647,6 +664,10 @@ def eval(cls, x, y):
         if isinstance(y, erf2inv) and y.args[0] == x:
             return y.args[1]
 
+        if x.is_zero or y.is_zero or x.is_extended_real and x.is_infinite or \
+                y.is_extended_real and y.is_infinite:
+            return erf(y) - erf(x)
+
         #Try to pull out -1 factor
         sign_x = x.could_extract_minus_sign()
         sign_y = y.could_extract_minus_sign()
@@ -770,6 +791,9 @@ def eval(cls, z):
         if isinstance(z, erf) and z.args[0].is_extended_real:
             return z.args[0]
 
+        if z.is_zero:
+            return S.Zero
+
         # Try to pull out factors of -1
         nz = z.extract_multiplicatively(-1)
         if nz is not None and (isinstance(nz, erf) and (nz.args[0]).is_extended_real):
@@ -778,6 +802,10 @@ def eval(cls, z):
     def _eval_rewrite_as_erfcinv(self, z, **kwargs):
        return erfcinv(1-z)
 
+    def _eval_is_zero(self):
+        if self.args[0].is_zero:
+            return True
+
 
 class erfcinv (Function):
     r"""
@@ -846,6 +874,9 @@ def eval(cls, z):
         elif z == 2:
             return S.NegativeInfinity
 
+        if z.is_zero:
+            return S.Infinity
+
     def _eval_rewrite_as_erfinv(self, z, **kwargs):
         return erfinv(1-z)
 
@@ -929,6 +960,18 @@ def eval(cls, x, y):
         elif y is S.Infinity:
             return erfinv(x)
 
+        if x.is_zero:
+            if y.is_zero:
+                return S.Zero
+            else:
+                return erfinv(y)
+        if y.is_zero:
+            return x
+
+    def _eval_is_zero(self):
+        x, y = self.args
+        if x.is_zero and y.is_zero:
+            return True
 
 ###############################################################################
 #################### EXPONENTIAL INTEGRALS ####################################
@@ -1028,6 +1071,9 @@ def eval(cls, z):
         elif z is S.NegativeInfinity:
             return S.Zero
 
+        if z.is_zero:
+            return S.NegativeInfinity
+
         nz, n = z.extract_branch_factor()
         if n:
             return Ei(nz) + 2*I*pi*n
@@ -1194,7 +1240,7 @@ def eval(cls, nu, z):
         # Extract branching information. This can be deduced from what is
         # explained in lowergamma.eval().
         z, n = z.extract_branch_factor()
-        if n == 0:
+        if n is S.Zero:
             return
         if nu.is_integer:
             if not nu > 0:
@@ -1377,6 +1423,8 @@ def eval(cls, z):
             return S.NegativeInfinity
         elif z is S.Infinity:
             return S.Infinity
+        if z.is_zero:
+            return S.Zero
 
     def fdiff(self, argindex=1):
         arg = self.args[0]
@@ -1424,6 +1472,10 @@ def _eval_rewrite_as_meijerg(self, z, **kwargs):
     def _eval_rewrite_as_tractable(self, z, **kwargs):
         return z * _eis(log(z))
 
+    def _eval_is_zero(self):
+        z = self.args[0]
+        if z.is_zero:
+            return True
 
 class Li(Function):
     r"""
@@ -1519,13 +1571,16 @@ class TrigonometricIntegral(Function):
 
     @classmethod
     def eval(cls, z):
-        if z == 0:
+        if z is S.Zero:
             return cls._atzero
         elif z is S.Infinity:
             return cls._atinf()
         elif z is S.NegativeInfinity:
             return cls._atneginf()
 
+        if z.is_zero:
+            return cls._atzero
+
         nz = z.extract_multiplicatively(polar_lift(I))
         if nz is None and cls._trigfunc(0) == 0:
             nz = z.extract_multiplicatively(I)
@@ -1672,10 +1727,16 @@ def _eval_rewrite_as_sinc(self, z, **kwargs):
         t = Symbol('t', Dummy=True)
         return Integral(sinc(t), (t, 0, z))
 
+    def _eval_is_zero(self):
+        z = self.args[0]
+        if z.is_zero:
+            return True
+
     def _sage_(self):
         import sage.all as sage
         return sage.sin_integral(self.args[0]._sage_())
 
+
 class Ci(TrigonometricIntegral):
     r"""
     Cosine integral.
@@ -1863,6 +1924,11 @@ def _eval_rewrite_as_expint(self, z, **kwargs):
         # XXX should we polarify z?
         return (E1(z) - E1(exp_polar(I*pi)*z))/2 - I*pi/2
 
+    def _eval_is_zero(self):
+        z = self.args[0]
+        if z.is_zero:
+            return True
+
     def _sage_(self):
         import sage.all as sage
         return sage.sinh_integral(self.args[0]._sage_())
@@ -1978,6 +2044,11 @@ class FresnelIntegral(Function):
 
     @classmethod
     def eval(cls, z):
+        # Values at positive infinities signs
+        # if any were extracted automatically
+        if z is S.Infinity:
+            return S.Half
+
         # Value at zero
         if z.is_zero:
             return S.Zero
@@ -2002,11 +2073,6 @@ def eval(cls, z):
         if changed:
             return prefact*cls(newarg)
 
-        # Values at positive infinities signs
-        # if any were extracted automatically
-        if z is S.Infinity:
-            return S.Half
-
     def fdiff(self, argindex=1):
         if argindex == 1:
             return self._trigfunc(S.Half*pi*self.args[0]**2)
@@ -2018,6 +2084,11 @@ def _eval_is_extended_real(self):
 
     _eval_is_finite = _eval_is_extended_real
 
+    def _eval_is_zero(self):
+        z = self.args[0]
+        if z.is_zero:
+            return True
+
     def _eval_conjugate(self):
         return self.func(self.args[0].conjugate())
 
@@ -2310,7 +2381,6 @@ class _erfs(Function):
     tractable for the Gruntz algorithm.
     """
 
-
     def _eval_aseries(self, n, args0, x, logx):
         from sympy import Order
         point = args0[0]
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 11608722f0..7e7ff701ac 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -303,7 +303,7 @@ def eval(cls, a, x):
         #    of both s and x), i.e.
         #    lowergamma(s, exp(2*I*pi*n)*x) = exp(2*pi*I*n*a)*lowergamma(a, x)
         from sympy import unpolarify, I
-        if x == 0:
+        if x is S.Zero:
             return S.Zero
         nx, n = x.extract_branch_factor()
         if a.is_integer and a.is_positive:
@@ -333,6 +333,9 @@ def eval(cls, a, x):
                 if not a.is_Integer:
                     return (-1)**(S.Half - a)*pi*erf(sqrt(x))/gamma(1 - a) + exp(-x)*Add(*[x**(k + a - 1)*gamma(a)/gamma(a + k) for k in range(1, Rational(3, 2) - a)])
 
+        if x.is_zero:
+            return S.Zero
+
     def _eval_evalf(self, prec):
         from mpmath import mp, workprec
         from sympy import Expr
@@ -346,9 +349,9 @@ def _eval_evalf(self, prec):
             return self
 
     def _eval_conjugate(self):
-        z = self.args[1]
-        if not z in (S.Zero, S.NegativeInfinity):
-            return self.func(self.args[0].conjugate(), z.conjugate())
+        x = self.args[1]
+        if x not in (S.Zero, S.NegativeInfinity):
+            return self.func(self.args[0].conjugate(), x.conjugate())
 
     def _eval_rewrite_as_uppergamma(self, s, x, **kwargs):
         return gamma(s) - uppergamma(s, x)
@@ -359,6 +362,11 @@ def _eval_rewrite_as_expint(self, s, x, **kwargs):
             return self
         return self.rewrite(uppergamma).rewrite(expint)
 
+    def _eval_is_zero(self):
+        x = self.args[1]
+        if x.is_zero:
+            return True
+
 
 class uppergamma(Function):
     r"""
@@ -488,6 +496,12 @@ def eval(cls, a, z):
                 if not a.is_Integer:
                     return (-1)**(S.Half - a) * pi*erfc(sqrt(z))/gamma(1-a) - z**a * exp(-z) * Add(*[z**k * gamma(a) / gamma(a+k+1) for k in range(S.Half - a)])
 
+        if a.is_zero and z.is_positive:
+            return -Ei(-z)
+
+        if z.is_zero and re(a).is_positive:
+            return gamma(a)
+
     def _eval_conjugate(self):
         z = self.args[1]
         if not z in (S.Zero, S.NegativeInfinity):
@@ -682,7 +696,7 @@ def eval(cls, n, z):
                 if z != nz:
                     return polygamma(n, nz)
 
-            if n == -1:
+            if n is S.NegativeOne:
                 return loggamma(z)
             else:
                 if z.is_Number:
@@ -694,6 +708,8 @@ def eval(cls, n, z):
                                 return S.Infinity
                             else:
                                 return S.Zero
+                        if n.is_zero:
+                            return S.Infinity
                     elif z.is_Integer:
                         if z.is_nonpositive:
                             return S.ComplexInfinity
@@ -703,7 +719,7 @@ def eval(cls, n, z):
                             elif n.is_odd:
                                 return (-1)**(n + 1)*factorial(n)*zeta(n + 1, z)
 
-        if n == 0:
+        if n.is_zero:
             if z is S.NaN:
                 return S.NaN
             elif z.is_Rational:
@@ -712,7 +728,7 @@ def eval(cls, n, z):
 
                 # only expand for small denominators to avoid creating long expressions
                 if q <= 5:
-                    return expand_func(polygamma(n, z, evaluate=False))
+                    return expand_func(polygamma(S.Zero, z, evaluate=False))
 
             elif z in (S.Infinity, S.NegativeInfinity):
                 return S.Infinity
diff --git a/sympy/functions/special/hyper.py b/sympy/functions/special/hyper.py
index 8aee18aa07..8162168e34 100644
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -178,9 +178,9 @@ class hyper(TupleParametersBase):
     """
 
 
-    def __new__(cls, ap, bq, z):
+    def __new__(cls, ap, bq, z, **kwargs):
         # TODO should we check convergence conditions?
-        return Function.__new__(cls, _prep_tuple(ap), _prep_tuple(bq), z)
+        return Function.__new__(cls, _prep_tuple(ap), _prep_tuple(bq), z, **kwargs)
 
     @classmethod
     def eval(cls, ap, bq, z):
@@ -442,7 +442,7 @@ class meijerg(TupleParametersBase):
     """
 
 
-    def __new__(cls, *args):
+    def __new__(cls, *args, **kwargs):
         if len(args) == 5:
             args = [(args[0], args[1]), (args[2], args[3]), args[4]]
         if len(args) != 3:
@@ -463,7 +463,7 @@ def tr(p):
                          "any b1, ..., bm by a positive integer")
 
         # TODO should we check convergence conditions?
-        return Function.__new__(cls, arg0, arg1, args[2])
+        return Function.__new__(cls, arg0, arg1, args[2], **kwargs)
 
     def fdiff(self, argindex=3):
         if argindex != 3:
diff --git a/sympy/functions/special/zeta_functions.py b/sympy/functions/special/zeta_functions.py
index 5eea2c4488..f82d3603c1 100644
--- a/sympy/functions/special/zeta_functions.py
+++ b/sympy/functions/special/zeta_functions.py
@@ -272,11 +272,11 @@ class polylog(Function):
     @classmethod
     def eval(cls, s, z):
         s, z = sympify((s, z))
-        if z == 1:
+        if z is S.One:
             return zeta(s)
-        elif z == -1:
+        elif z is S.NegativeOne:
             return -dirichlet_eta(s)
-        elif z == 0:
+        elif z is S.Zero:
             return S.Zero
         elif s == 2:
             if z == S.Half:
@@ -291,17 +291,32 @@ def eval(cls, s, z):
                 return pi**2/15 - log((sqrt(5)-1)/2)**2
             elif z == (sqrt(5) - 1)/2:
                 return pi**2/10 - log((sqrt(5)-1)/2)**2
-        # For s = 0 or -1 use explicit formulas to evaluate, but
-        # automatically expanding polylog(1, z) to -log(1-z) seems undesirable
-        # for summation methods based on hypergeometric functions
-        elif s == 0:
-            return z/(1 - z)
-        elif s == -1:
-            return z/(1 - z)**2
+
+        if z.is_zero:
+            return S.Zero
+
+        # Make an effort to determine if z is 1 to avoid replacing into
+        # expression with singularity
+        zone = z.equals(S.One)
+
+        if zone:
+            return zeta(s)
+        elif zone is False:
+            # For s = 0 or -1 use explicit formulas to evaluate, but
+            # automatically expanding polylog(1, z) to -log(1-z) seems
+            # undesirable for summation methods based on hypergeometric
+            # functions
+            if s is S.Zero:
+                return z/(1 - z)
+            elif s is S.NegativeOne:
+                return z/(1 - z)**2
+            if s.is_zero:
+                return z/(1 - z)
+
         # polylog is branched, but not over the unit disk
         from sympy.functions.elementary.complexes import (Abs, unpolarify,
-            polar_lift)
-        if z.has(exp_polar, polar_lift) and (Abs(z) <= S.One) == True:
+                                                          polar_lift)
+        if z.has(exp_polar, polar_lift) and (zone or (Abs(z) <= S.One) == True):
             return cls(s, unpolarify(z))
 
     def fdiff(self, argindex=1):
@@ -326,6 +341,11 @@ def _eval_expand_func(self, **hints):
             return expand_mul(start).subs(u, z)
         return polylog(s, z)
 
+    def _eval_is_zero(self):
+        z = self.args[1]
+        if z.is_zero:
+            return True
+
 ###############################################################################
 ###################### HURWITZ GENERALIZED ZETA FUNCTION ######################
 ###############################################################################
@@ -476,7 +496,8 @@ def eval(cls, z, a_=None):
                     return zeta + harmonic(abs(a), z)
                 else:
                     return zeta - harmonic(a - 1, z)
-
+        if z.is_zero:
+            return S.Half - a
 
     def _eval_rewrite_as_dirichlet_eta(self, s, a=1, **kwargs):
         if a != 1:
@@ -601,5 +622,14 @@ def eval(cls, n, a=None):
                 return S.ComplexInfinity
             elif not n.is_Integer:
                 return S.ComplexInfinity
-            elif n == 0 and a in [None, 1]:
+            elif n is S.Zero and a in [None, 1]:
                 return S.EulerGamma
+
+        if n.is_extended_negative:
+            return S.ComplexInfinity
+
+        if n.is_zero and a in [None, 1]:
+            return S.EulerGamma
+
+        if n.is_integer == False:
+            return S.ComplexInfinity
diff --git a/sympy/stats/frv_types.py b/sympy/stats/frv_types.py
index c6d357e29b..1e01e0334a 100644
--- a/sympy/stats/frv_types.py
+++ b/sympy/stats/frv_types.py
@@ -433,7 +433,7 @@ def BetaBinomial(name, n, alpha, beta):
 
     >>> X = BetaBinomial('X', 2, 1, 1)
     >>> density(X).dict
-    {0: beta(1, 3)/beta(1, 1), 1: 2*beta(2, 2)/beta(1, 1), 2: beta(3, 1)/beta(1, 1)}
+    {0: 1/3, 1: 2*beta(2, 2), 2: 1/3}
 
     References
     ==========


Test patch:
diff --git a/sympy/functions/elementary/tests/test_trigonometric.py b/sympy/functions/elementary/tests/test_trigonometric.py
index f7012653b3..bba1278acc 100644
--- a/sympy/functions/elementary/tests/test_trigonometric.py
+++ b/sympy/functions/elementary/tests/test_trigonometric.py
@@ -3,7 +3,7 @@
         cosh, atan2, exp, log, asinh, acoth, atanh, O, cancel, Matrix, re, im,
         Float, Pow, gcd, sec, csc, cot, diff, simplify, Heaviside, arg,
         conjugate, series, FiniteSet, asec, acsc, Mul, sinc, jn,
-        AccumBounds, Interval, ImageSet, Lambda, besselj)
+        AccumBounds, Interval, ImageSet, Lambda, besselj, Add)
 from sympy.core.compatibility import range
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
@@ -147,6 +147,11 @@ def test_sin():
             e = abs( float(sin(x)) - sin(float(x)) )
             assert e < 1e-12
 
+    assert sin(0, evaluate=False).is_zero is True
+    assert sin(k*pi, evaluate=False).is_zero is None
+
+    assert sin(Add(1, -1, evaluate=False), evaluate=False).is_zero is True
+
 
 def test_sin_cos():
     for d in [1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 24, 30, 40, 60, 120]:  # list is not exhaustive...
diff --git a/sympy/functions/special/tests/test_hyper.py b/sympy/functions/special/tests/test_hyper.py
index 69281d61b0..7afcf819b6 100644
--- a/sympy/functions/special/tests/test_hyper.py
+++ b/sympy/functions/special/tests/test_hyper.py
@@ -1,7 +1,7 @@
 from sympy import (hyper, meijerg, S, Tuple, pi, I, exp, log, Rational,
                    cos, sqrt, symbols, oo, Derivative, gamma, O, appellf1)
-from sympy.series.limits import limit
 from sympy.abc import x, z, k
+from sympy.series.limits import limit
 from sympy.utilities.pytest import raises, slow
 from sympy.utilities.randtest import (
     random_complex_number as randcplx,
@@ -47,6 +47,10 @@ def test_hyper():
     assert hyper([polar_lift(z)], [polar_lift(k)], polar_lift(x)) == \
         hyper([z], [k], polar_lift(x))
 
+    # hyper does not automatically evaluate anyway, but the test is to make
+    # sure that the evaluate keyword is accepted
+    assert hyper((1, 2), (1,), z, evaluate=False).func is hyper
+
 
 def test_expand_func():
     # evaluation at 1 of Gauss' hypergeometric function:
@@ -350,6 +354,11 @@ def test_appellf1():
     assert appellf1(a, b1, b1, c, y, x) == appellf1(a, b1, b1, c, x, y)
     assert appellf1(a, b1, b2, c, S.Zero, S.Zero) is S.One
 
+    f = appellf1(a, b1, b2, c, S.Zero, S.Zero, evaluate=False)
+    assert f.func is appellf1
+    assert f.doit() is S.One
+
+
 def test_derivative_appellf1():
     from sympy import diff
     a, b1, b2, c, x, y, z = symbols('a b1 b2 c x y z')
diff --git a/sympy/functions/special/tests/test_zeta_functions.py b/sympy/functions/special/tests/test_zeta_functions.py
index 3eaf4d2c81..573d6c59bd 100644
--- a/sympy/functions/special/tests/test_zeta_functions.py
+++ b/sympy/functions/special/tests/test_zeta_functions.py
@@ -167,6 +167,9 @@ def test_polylog_values():
             assert tn(polylog(s, z), polylog(s, z, evaluate=False), z,
                 a=2, b=-2, c=5, d=2)
 
+    from sympy import Integral
+    assert polylog(0, Integral(1, (x, 0, 1))) == -S.Half
+
 
 def test_lerchphi_expansion():
     assert myexpand(lerchphi(1, s, a), zeta(s, a))