Repository: sympy/sympy

Issue:
`harmonic(n,m).evalf()` fails on non-integer `m`
```
>>> from sympy import *
>>> harmonic(1,pi).evalf()
0
>>> harmonic(2,pi).evalf()
0
```

But I thought they would be $1$ and $1+2^{-\pi}$?

Discussion:

Patch:
diff --git a/sympy/functions/combinatorial/numbers.py b/sympy/functions/combinatorial/numbers.py
index eda3546f56..a4cbb74cdb 100644
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -14,7 +14,7 @@
 from sympy.core.cache import cacheit
 from sympy.core.evalf import pure_complex
 from sympy.core.expr import Expr
-from sympy.core.function import Function, expand_mul
+from sympy.core.function import ArgumentIndexError, Function, expand_mul
 from sympy.core.logic import fuzzy_not
 from sympy.core.mul import Mul
 from sympy.core.numbers import E, pi, oo, Rational, Integer
@@ -781,6 +781,12 @@ class harmonic(Function):
       of order `m`, `\operatorname{H}_{n,m}`, where
       ``harmonic(n) == harmonic(n, 1)``
 
+    This function can be extended to complex `n` and `m` where `n` is not a
+    negative integer as
+
+    .. math:: \operatorname{H}_{n,m} = \begin{cases} \zeta(m) - \zeta(m, n+1)
+            & m \ne 1 \\ \psi(n+1) + \gamma & m = 1 \end{cases}
+
     Examples
     ========
 
@@ -829,7 +835,7 @@ class harmonic(Function):
     We can rewrite harmonic numbers in terms of polygamma functions:
 
     >>> from sympy import digamma, polygamma
-    >>> m = Symbol("m")
+    >>> m = Symbol("m", integer=True, positive=True)
 
     >>> harmonic(n).rewrite(digamma)
     polygamma(0, n + 1) + EulerGamma
@@ -841,7 +847,7 @@ class harmonic(Function):
     polygamma(2, n + 1)/2 - polygamma(2, 1)/2
 
     >>> harmonic(n,m).rewrite(polygamma)
-    (-1)**m*(polygamma(m - 1, 1) - polygamma(m - 1, n + 1))/factorial(m - 1)
+    (-1)**m*(polygamma(m - 1, 1) - polygamma(m - 1, n + 1))/gamma(m)
 
     Integer offsets in the argument can be pulled out:
 
@@ -866,12 +872,8 @@ class harmonic(Function):
     >>> limit(harmonic(n, 3), n, oo)
     -polygamma(2, 1)/2
 
-    However we cannot compute the general relation yet:
-
-    >>> limit(harmonic(n, m), n, oo)
-    harmonic(oo, m)
-
-    which equals ``zeta(m)`` for ``m > 1``.
+    >>> limit(harmonic(n, m+1), n, oo) # does not hold for m == 1
+    zeta(m + 1)
 
     See Also
     ========
@@ -898,34 +900,31 @@ def eval(cls, n, m=None):
             return cls(n)
         if m is None:
             m = S.One
-
-        if m.is_zero:
+        if n.is_zero:
+            return S.Zero
+        elif m.is_zero:
             return n
-
-        if n is S.Infinity:
+        elif n is S.Infinity:
             if m.is_negative:
                 return S.NaN
             elif is_le(m, S.One):
                 return S.Infinity
             elif is_gt(m, S.One):
                 return zeta(m)
-            else:
-                return
-
-        if n == 0:
-            return S.Zero
-
-        if n.is_Integer and n.is_nonnegative:
-            if m not in cls._functions:
+        elif n.is_Integer:
+            if n.is_negative:
+                return S.NaN
+            elif m not in cls._functions:
                 @recurrence_memo([0])
                 def f(n, prev):
                     return prev[-1] + S.One / n**m
                 cls._functions[m] = f
             return cls._functions[m](int(n))
 
-    def _eval_rewrite_as_polygamma(self, n, m=1, **kwargs):
-        from sympy.functions.special.gamma_functions import polygamma
-        return S.NegativeOne**m/factorial(m - 1) * (polygamma(m - 1, 1) - polygamma(m - 1, n + 1))
+    def _eval_rewrite_as_polygamma(self, n, m=S.One, **kwargs):
+        from sympy.functions.special.gamma_functions import gamma, polygamma
+        if m.is_integer and m.is_positive:
+            return S.NegativeOne**m * (polygamma(m-1, 1) - polygamma(m-1, n+1)) / gamma(m)
 
     def _eval_rewrite_as_digamma(self, n, m=1, **kwargs):
         from sympy.functions.special.gamma_functions import polygamma
@@ -942,6 +941,12 @@ def _eval_rewrite_as_Sum(self, n, m=None, **kwargs):
             m = S.One
         return Sum(k**(-m), (k, 1, n))
 
+    def _eval_rewrite_as_zeta(self, n, m=S.One, **kwargs):
+        from sympy.functions.special.zeta_functions import zeta
+        from sympy.functions.special.gamma_functions import digamma
+        return Piecewise((digamma(n + 1) + S.EulerGamma, Eq(m, 1)),
+                         (zeta(m) - zeta(m, n+1), True))
+
     def _eval_expand_func(self, **hints):
         from sympy.concrete.summations import Sum
         n = self.args[0]
@@ -979,13 +984,37 @@ def _eval_expand_func(self, **hints):
         return self
 
     def _eval_rewrite_as_tractable(self, n, m=1, limitvar=None, **kwargs):
+        from sympy.functions.special.zeta_functions import zeta
         from sympy.functions.special.gamma_functions import polygamma
-        return self.rewrite(polygamma).rewrite("tractable", deep=True)
+        pg = self.rewrite(polygamma)
+        if not isinstance(pg, harmonic):
+            return pg.rewrite("tractable", deep=True)
+        return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
 
     def _eval_evalf(self, prec):
-        from sympy.functions.special.gamma_functions import polygamma
-        if all(i.is_number for i in self.args):
-            return self.rewrite(polygamma)._eval_evalf(prec)
+        if not all(x.is_number for x in self.args):
+            return
+        n = self.args[0]._to_mpmath(prec)
+        m = (self.args[1] if len(self.args) > 1 else S.One)._to_mpmath(prec)
+        if mp.isint(n) and n < 0:
+            return S.NaN
+        with workprec(prec):
+            if m == 1:
+                res = mp.harmonic(n)
+            else:
+                res = mp.zeta(m) - mp.zeta(m, n+1)
+        return Expr._from_mpmath(res, prec)
+
+    def fdiff(self, argindex=1):
+        from sympy.functions.special.zeta_functions import zeta
+        if len(self.args) == 2:
+            n, m = self.args
+        else:
+            n, m = self.args + (1,)
+        if argindex == 1:
+            return m * zeta(m+1, n+1)
+        else:
+            raise ArgumentIndexError
 
 
 #----------------------------------------------------------------------------#
diff --git a/sympy/functions/special/zeta_functions.py b/sympy/functions/special/zeta_functions.py
index 71fdb6a34e..ce067a8c86 100644
--- a/sympy/functions/special/zeta_functions.py
+++ b/sympy/functions/special/zeta_functions.py
@@ -506,6 +506,8 @@ def eval(cls, s, a=None):
             return S.ComplexInfinity
         elif s is S.Infinity:
             return S.One
+        elif a is S.Infinity:
+            return S.Zero
 
         sint = s.is_Integer
         if a is None:
diff --git a/sympy/printing/numpy.py b/sympy/printing/numpy.py
index 5f6186e411..d9544a4006 100644
--- a/sympy/printing/numpy.py
+++ b/sympy/printing/numpy.py
@@ -404,6 +404,9 @@ def _print_bernoulli(self, expr):
         # scipy's bernoulli is inconsistent with SymPy's so rewrite
         return self._print(expr._eval_rewrite_as_zeta(*expr.args))
 
+    def _print_harmonic(self, expr):
+        return self._print(expr._eval_rewrite_as_zeta(*expr.args))
+
     def _print_Integral(self, e):
         integration_vars, limits = _unpack_integral_limits(e)
 


Test patch:
diff --git a/sympy/functions/combinatorial/tests/test_comb_numbers.py b/sympy/functions/combinatorial/tests/test_comb_numbers.py
index 376e4a1761..91d9091e94 100644
--- a/sympy/functions/combinatorial/tests/test_comb_numbers.py
+++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py
@@ -12,6 +12,7 @@
 from sympy.functions.elementary.integers import floor
 from sympy.polys.polytools import cancel
 from sympy.series.limits import limit
+from sympy.series.order import O
 from sympy.functions import (
     bernoulli, harmonic, bell, fibonacci, tribonacci, lucas, euler, catalan,
     genocchi, partition, motzkin, binomial, gamma, sqrt, cbrt, hyper, log, digamma,
@@ -21,8 +22,7 @@
 from sympy.core.expr import unchanged
 from sympy.core.numbers import GoldenRatio, Integer
 
-from sympy.testing.pytest import (XFAIL, raises, nocache_fail,
-                                  warns_deprecated_sympy)
+from sympy.testing.pytest import raises, nocache_fail, warns_deprecated_sympy
 from sympy.abc import x
 
 
@@ -246,6 +246,8 @@ def test_harmonic():
     assert harmonic(oo, 1 + ip) is zeta(1 + ip)
 
     assert harmonic(0, m) == 0
+    assert harmonic(-1, 2) is S.NaN
+    assert harmonic(-2, n) is S.NaN
 
 
 def test_harmonic_rational():
@@ -315,35 +317,47 @@ def test_harmonic_rational():
 def test_harmonic_evalf():
     assert str(harmonic(1.5).evalf(n=10)) == '1.280372306'
     assert str(harmonic(1.5, 2).evalf(n=10)) == '1.154576311'  # issue 7443
+    assert str(harmonic(4.0, -3).evalf(n=10)) == '100.0000000'
+    assert str(harmonic(7.0, 1.0).evalf(n=10)) == '2.592857143'
+    assert str(harmonic(1, pi).evalf(n=10)) == '1.000000000'
+    assert str(harmonic(2, pi).evalf(n=10)) == '1.113314732'
+    assert str(harmonic(1000.0, pi).evalf(n=10)) == '1.176241563'
+    assert str(harmonic(I).evalf(n=10)) == '0.6718659855 + 1.076674047*I'
+    assert str(harmonic(I, I).evalf(n=10)) == '-0.3970915266 + 1.9629689*I'
 
+    assert harmonic(-1.0, 1).evalf() is S.NaN
+    assert harmonic(-2.0, 2.0).evalf() is S.NaN
 
 def test_harmonic_rewrite():
     n = Symbol("n")
-    m = Symbol("m")
+    m = Symbol("m", integer=True, positive=True)
 
     assert harmonic(n).rewrite(digamma) == polygamma(0, n + 1) + EulerGamma
-    assert harmonic(n).rewrite(trigamma) ==  polygamma(0, n + 1) + EulerGamma
-    assert harmonic(n).rewrite(polygamma) ==  polygamma(0, n + 1) + EulerGamma
+    assert harmonic(n).rewrite(trigamma) == polygamma(0, n + 1) + EulerGamma
+    assert harmonic(n).rewrite(polygamma) == polygamma(0, n + 1) + EulerGamma
 
     assert harmonic(n,3).rewrite(polygamma) == polygamma(2, n + 1)/2 - polygamma(2, 1)/2
-    assert harmonic(n,m).rewrite(polygamma) == (-1)**m*(polygamma(m - 1, 1) - polygamma(m - 1, n + 1))/factorial(m - 1)
+    assert harmonic(n,m).rewrite(polygamma) == (-1)**m * (polygamma(m-1, 1) - polygamma(m-1, n+1)) / gamma(m)
 
     assert expand_func(harmonic(n+4)) == harmonic(n) + 1/(n + 4) + 1/(n + 3) + 1/(n + 2) + 1/(n + 1)
     assert expand_func(harmonic(n-4)) == harmonic(n) - 1/(n - 1) - 1/(n - 2) - 1/(n - 3) - 1/n
 
     assert harmonic(n, m).rewrite("tractable") == harmonic(n, m).rewrite(polygamma)
+    assert harmonic(n, x).rewrite("tractable") == zeta(x) - zeta(x, n+1)
 
     _k = Dummy("k")
     assert harmonic(n).rewrite(Sum).dummy_eq(Sum(1/_k, (_k, 1, n)))
     assert harmonic(n, m).rewrite(Sum).dummy_eq(Sum(_k**(-m), (_k, 1, n)))
 
 
-@XFAIL
-def test_harmonic_limit_fail():
-    n = Symbol("n")
-    m = Symbol("m")
-    # For m > 1:
-    assert limit(harmonic(n, m), n, oo) == zeta(m)
+def test_harmonic_calculus():
+    y = Symbol("y", positive=True)
+    assert harmonic(x, 1).limit(x, 0) == 0
+    assert harmonic(x, y).limit(x, 0) == 0
+    assert harmonic(x, 1).series(x, y, 2) == \
+            harmonic(y) + (x - y)*zeta(2, y + 1) + O((x - y)**2, (x, y))
+    assert limit(harmonic(x, y), x, oo) == harmonic(oo, y)
+    assert limit(harmonic(x, y+1), x, oo) == zeta(y+1)
 
 
 def test_euler():
diff --git a/sympy/utilities/tests/test_lambdify.py b/sympy/utilities/tests/test_lambdify.py
index f06c210c01..80a7dcb286 100644
--- a/sympy/utilities/tests/test_lambdify.py
+++ b/sympy/utilities/tests/test_lambdify.py
@@ -11,7 +11,7 @@
 from sympy.core.singleton import S
 from sympy.core.symbol import (Dummy, symbols)
 from sympy.functions.combinatorial.factorials import (RisingFactorial, factorial)
-from sympy.functions.combinatorial.numbers import bernoulli
+from sympy.functions.combinatorial.numbers import bernoulli, harmonic
 from sympy.functions.elementary.complexes import Abs
 from sympy.functions.elementary.exponential import exp, log
 from sympy.functions.elementary.hyperbolic import acosh
@@ -1476,6 +1476,16 @@ def test_scipy_bernoulli():
     assert bern(1) == 0.5
 
 
+def test_scipy_harmonic():
+    if not scipy:
+        skip("scipy not installed")
+
+    hn = lambdify((x,), harmonic(x), modules='scipy')
+    assert hn(2) == 1.5
+    hnm = lambdify((x, y), harmonic(x, y), modules='scipy')
+    assert hnm(2, 2) == 1.25
+
+
 def test_cupy_array_arg():
     if not cupy:
         skip("CuPy not installed")