Repository: sympy/sympy
Pull request: polygamma improvements (evaluation, real and positivity determination)
Closes #17350. Fixes #12569.

#### Brief description of what is fixed or changed

- Add `_eval_evalf` to prevent evaluation when the first argument is not a nonnegative integer.
- Remove `_eval_is_real` since the logic was wrong.
- Update `_eval_is_positive` and `_eval_is_negative` to avoid doing comparisons with `>`.

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->


Problem statement:
problem with polygamma or im
N(polygamma(0,I))
yields a correct result 
0.094650320622477 + 2.07667404746858*I
but
im(polygamma(0,I))
yields zero.

Should polygamma(a, b) raise NotImplementedError for noninteger (numeric) a?
mpmath's implementation of polygamma is valid only for integer first arguments. SymPy's implementation does not care, but this leads to many issues if you actually try to do things with such objects:

```
>>> (I*polygamma(I, pi)).as_real_imag()
Traceback (most recent call last):
  File "/home/eward/se2/sympy/core/assumptions.py", line 262, in getit
    return self._assumptions[fact]
KeyError: 'zero'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/eward/se2/sympy/core/assumptions.py", line 262, in getit
    return self._assumptions[fact]
KeyError: 'finite'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/eward/se2/sympy/core/assumptions.py", line 262, in getit
    return self._assumptions[fact]
KeyError: 'extended_positive'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/eward/se2/sympy/core/evalf.py", line 1308, in evalf
    rf = evalf_table[x.func]
KeyError: polygamma

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/eward/se2/sympy/core/mul.py", line 798, in as_real_imag
    if i.is_zero:
  File "/home/eward/se2/sympy/core/assumptions.py", line 266, in getit
    return _ask(fact, self)
  File "/home/eward/se2/sympy/core/assumptions.py", line 321, in _ask
    _ask(pk, obj)
  File "/home/eward/se2/sympy/core/assumptions.py", line 309, in _ask
    a = evaluate(obj)
  File "/home/eward/se2/sympy/core/expr.py", line 864, in _eval_is_negative
    finite = self.is_finite
  File "/home/eward/se2/sympy/core/assumptions.py", line 266, in getit
    return _ask(fact, self)
  File "/home/eward/se2/sympy/core/assumptions.py", line 321, in _ask
    _ask(pk, obj)
  File "/home/eward/se2/sympy/core/assumptions.py", line 309, in _ask
    a = evaluate(obj)
  File "/home/eward/se2/sympy/core/expr.py", line 857, in _eval_is_positive
    extended_positive = self.is_extended_positive
  File "/home/eward/se2/sympy/core/assumptions.py", line 266, in getit
    return _ask(fact, self)
  File "/home/eward/se2/sympy/core/assumptions.py", line 309, in _ask
    a = evaluate(obj)
  File "/home/eward/se2/sympy/core/expr.py", line 882, in _eval_is_extended_positive
    n2 = self._eval_evalf(2)
  File "/home/eward/se2/sympy/core/function.py", line 565, in _eval_evalf
    args = [arg._to_mpmath(prec + 5) for arg in self.args]
  File "/home/eward/se2/sympy/core/function.py", line 565, in <listcomp>
    args = [arg._to_mpmath(prec + 5) for arg in self.args]
  File "/home/eward/se2/sympy/core/evalf.py", line 1489, in _to_mpmath
    re, im, _, _ = evalf(self, prec, {})
  File "/home/eward/se2/sympy/core/evalf.py", line 1314, in evalf
    xe = x._eval_evalf(prec)
  File "/home/eward/se2/sympy/core/function.py", line 590, in _eval_evalf
    v = func(*args)
  File "/usr/local/lib/python3.6/dist-packages/mpmath/ctx_mp.py", line 265, in psi
    m = int(m)
TypeError: int() argument must be a string, a bytes-like object or a number, not 'mpc'
>>> (tanh(polygamma(I, 1))).rewrite(exp)
...
    m = int(m)
TypeError: int() argument must be a string, a bytes-like object or a number, not 'mpc'
>>> (I / polygamma(I, 4)).rewrite(exp)
...
    raise TypeError("Invalid comparison of complex %s" % me)
TypeError: Invalid comparison of complex I
```

Hints from issue discussion:
`polygamma` is real for real arguments. Hence the imaginary part is zero.
```
>>> N(polygamma(0, 1))
-0.577215664901533
>>> im(polygamma(0, 1))
0
```

Please note that the second argument is I (imaginary unit) not 1 (one).
I see, now that you point that out. (My eyesight is not what it used to be..)
It seems that the error results from this:
```
>>> polygamma(0, I).is_real
True
```
There a bug in `polygamma._eval_is_real`. It does not check the proper argument.
I think we can set `is_number` to False to indicate to things that they shouldn't try to evalf. 
There is a more generalized definition for the polygamma function 
http://dauns01.math.tulane.edu/~vhm/papers_html/genoff.pdf
Can it be used here?
Yes, but that would be on the mpmath side. If that's not changed then SymPy needs to be more careful about evaluation.

Commit message: polygamma improvements: prevent bad evaluation, update sign, real determination

Patch to review:
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index ee79a3597770..57d1cb3a0ef4 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -583,6 +583,12 @@ class polygamma(Function):
     .. [4] http://functions.wolfram.com/GammaBetaErf/PolyGamma2/
     """
 
+    def _eval_evalf(self, prec):
+        n = self.args[0]
+        # the mpmath polygamma implementation valid only for nonnegative integers
+        if n.is_real and n.is_number:
+            if (n.is_integer or n == int(n)) and n.is_nonnegative:
+                return super(polygamma, self)._eval_evalf(prec)
 
     def fdiff(self, argindex=2):
         if argindex == 2:
@@ -592,16 +598,13 @@ def fdiff(self, argindex=2):
             raise ArgumentIndexError(self, argindex)
 
     def _eval_is_positive(self):
-        if self.args[1].is_positive and (self.args[0] > 0) == True:
+        if self.args[0].is_positive and self.args[1].is_positive:
             return self.args[0].is_odd
 
     def _eval_is_negative(self):
-        if self.args[1].is_positive and (self.args[0] > 0) == True:
+        if self.args[0].is_positive and self.args[1].is_positive:
             return self.args[0].is_even
 
-    def _eval_is_real(self):
-        return self.args[0].is_real
-
     def _eval_aseries(self, n, args0, x, logx):
         from sympy import Order
         if args0[1] != oo or not \
diff --git a/sympy/functions/special/tests/test_gamma_functions.py b/sympy/functions/special/tests/test_gamma_functions.py
index 1c94267bee41..3c6f1eb62eca 100644
--- a/sympy/functions/special/tests/test_gamma_functions.py
+++ b/sympy/functions/special/tests/test_gamma_functions.py
@@ -2,7 +2,7 @@
     Symbol, Dummy, gamma, I, oo, nan, zoo, factorial, sqrt, Rational,
     multigamma, log, polygamma, EulerGamma, pi, uppergamma, S, expand_func,
     loggamma, sin, cos, O, lowergamma, exp, erf, erfc, exp_polar, harmonic,
-    zeta, conjugate)
+    zeta, conjugate, im, re, tanh, Abs)
 
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
@@ -275,6 +275,38 @@ def t(m, n):
     # Test a bug
     assert polygamma(0, -x).expand(func=True) == polygamma(0, -x)
 
+    assert polygamma(2, 2.5).is_positive == False
+    assert polygamma(2, -2.5).is_positive == False
+    assert polygamma(3, 2.5).is_positive == True
+    assert polygamma(3, -2.5).is_positive is None
+    assert polygamma(-2, -2.5).is_positive is None
+    assert polygamma(-3, -2.5).is_positive is None
+
+    assert polygamma(2, 2.5).is_negative == True
+    assert polygamma(3, 2.5).is_negative == False
+    assert polygamma(3, -2.5).is_negative == False
+    assert polygamma(2, -2.5).is_negative is None
+    assert polygamma(-2, -2.5).is_negative is None
+    assert polygamma(-3, -2.5).is_negative is None
+
+    assert polygamma(I, 2).is_positive is None
+    assert polygamma(I, 3).is_negative is None
+
+    # issue 17350
+    assert polygamma(pi, 3).evalf() == polygamma(pi, 3)
+    assert (I*polygamma(I, pi)).as_real_imag() == \
+           (-im(polygamma(I, pi)), re(polygamma(I, pi)))
+    assert (tanh(polygamma(I, 1))).rewrite(exp) == \
+           (exp(polygamma(I, 1)) - exp(-polygamma(I, 1)))/(exp(polygamma(I, 1)) + exp(-polygamma(I, 1)))
+    assert (I / polygamma(I, 4)).rewrite(exp) == \
+           I*sqrt(re(polygamma(I, 4))**2 + im(polygamma(I, 4))**2)\
+           /((re(polygamma(I, 4)) + I*im(polygamma(I, 4)))*Abs(polygamma(I, 4)))
+    assert unchanged(polygamma, 2.3, 1.0)
+
+    # issue 12569
+    assert unchanged(im, polygamma(0, I))
+    assert polygamma(0, I).is_real is None
+
 
 def test_polygamma_expand_func():
     assert polygamma(0, x).expand(func=True) == polygamma(0, x)


Code review comments:
diff_hunk: @@ -583,6 +583,12 @@ class polygamma(Function):
     .. [4] http://functions.wolfram.com/GammaBetaErf/PolyGamma2/
     """
 
+    def _eval_evalf(self, prec):
+        n = self.args[0]
+        # the mpmath polygamma implementation valid only for nonnegative integers
+        if n.is_real and n.is_number:
path: sympy/functions/special/gamma_functions.py
text: @user1:
How about changing the order? I think that `n.is_number` will generally take less time than `n.is_real` (though I may be wrong).

@user2:
It should be correct that `n.is_number` is faster as it only has to check if there are any free symbols and not worry about the assumptions on the symbols.