diff --git a/src/sage/functions/all.py b/src/sage/functions/all.py
index 0eb6add9839..53b0cb654bf 100644
--- a/src/sage/functions/all.py
+++ b/src/sage/functions/all.py
@@ -26,7 +26,7 @@
                     arg, real_part, real, frac,
                     imag_part, imag, imaginary, conjugate)
 
-from .log import (exp, exp_polar, log, ln, polylog, dilog, lambert_w)
+from .log import (exp, exp_polar, log, ln, polylog, dilog, lambert_w, harmonic_number)
 
 
 from .transcendental import (zeta, zetaderiv, zeta_symmetric, hurwitz_zeta,
diff --git a/src/sage/functions/log.py b/src/sage/functions/log.py
index b92d87d379b..95a6ae6b9c2 100644
--- a/src/sage/functions/log.py
+++ b/src/sage/functions/log.py
@@ -18,6 +18,9 @@
 from sage.rings.real_double import RDF
 from sage.rings.complex_double import CDF
 from sage.rings.integer import Integer
+from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+from sage.rings.rational import Rational
 
 class Function_exp(GinacFunction):
     def __init__(self):
@@ -813,6 +816,7 @@ def _print_latex_(self, n, z):
 
 lambert_w = Function_lambert_w()
 
+
 class Function_exp_polar(BuiltinFunction):
     def __init__(self):
         r"""
@@ -835,7 +839,7 @@ def __init__(self):
         function::
 
             sage: exp_polar(pi*I/2)
-            I       
+            I
             sage: x = var('x', domain='real')
             sage: exp_polar(-1/2*I*pi + x)
             e^(-1/2*I*pi + x)
@@ -916,3 +920,355 @@ def _eval_(self, z):
             return None
 
 exp_polar = Function_exp_polar()
+
+
+class Function_harmonic_number_generalized(BuiltinFunction):
+    r"""
+    Harmonic and generalized harmonic number functions,
+    defined by:
+
+    .. math::
+
+        H_{n}=H_{n,1}=\sum_{k=1}^n\frac1k
+        
+        H_{n,m}=\sum_{k=1}^n\frac1{k^m}
+
+    They are also well-defined for complex argument, through:
+
+    .. math::
+
+        H_{s}=\int_0^1\frac{1-x^s}{1-x}
+
+        H_{s,m}=\zeta(m)-\zeta(m,s-1)
+
+    If called with a single argument, that argument is ``s`` and ``m`` is
+    assumed to be 1 (the normal harmonic numbers ``H_s``).
+
+    ALGORITHM:
+
+    Numerical evaluation is handled using the mpmath and FLINT libraries.
+
+    REFERENCES:
+
+    - :wikipedia:`/Harmonic_number`
+
+    EXAMPLES:
+
+    Evaluation of integer, rational, or complex argument::
+
+        sage: harmonic_number(5)
+        137/60
+        sage: harmonic_number(3,3)
+        251/216
+        sage: harmonic_number(5/2)
+        -2*log(2) + 46/15
+        sage: harmonic_number(3.,3)
+        zeta(3) - 0.0400198661225573
+        sage: harmonic_number(3.,3.)
+        1.16203703703704
+        sage: harmonic_number(3,3).n(200)
+        1.16203703703703703703703...
+        sage: harmonic_number(1+I,5)
+        harmonic_number(I + 1, 5)
+        sage: harmonic_number(5,1.+I)
+        1.57436810798989 - 1.06194728851357*I
+
+    Solutions to certain sums are returned in terms of harmonic numbers::
+
+        sage: k=var('k')
+        sage: sum(1/k^7,k,1,x)
+        harmonic_number(x, 7)
+
+    Check the defining integral at a random integer::
+
+        sage: n=randint(10,100)
+        sage: bool(SR(integrate((1-x^n)/(1-x),x,0,1)) == harmonic_number(n))
+        True
+
+    There are several special values which are automatically simplified::
+
+        sage: harmonic_number(0)
+        0
+        sage: harmonic_number(1)
+        1
+        sage: harmonic_number(x,1)
+        harmonic_number(x)
+
+    Arguments are swapped with respect to the same functions in
+    Maxima::
+
+        sage: maxima(harmonic_number(x,2)) # maxima expect interface
+        gen_harmonic_number(2,_SAGE_VAR_x)
+        sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefms
+        sage: sefms('gen_harmonic_number(3,x)')
+        harmonic_number(x, 3)
+        sage: from sage.interfaces.maxima_lib import maxima_lib, max_to_sr
+        sage: c=maxima_lib(harmonic_number(x,2)); c
+        gen_harmonic_number(2,_SAGE_VAR_x)
+        sage: max_to_sr(c.ecl())
+        harmonic_number(x, 2)
+    """
+
+    def __init__(self):
+        r"""
+        EXAMPLES::
+
+            sage: loads(dumps(harmonic_number(x,5)))
+            harmonic_number(x, 5)
+        """
+        BuiltinFunction.__init__(self, "harmonic_number", nargs=2)
+
+    def __call__(self, z, m=1, **kwds):
+        r"""
+        Custom call method allows the user to pass one argument or two. If
+        one argument is passed, we assume it is ``z`` and that ``m=1``.
+
+        EXAMPLES::
+
+            sage: harmonic_number(x)
+            harmonic_number(x)
+            sage: harmonic_number(x,1)
+            harmonic_number(x)
+            sage: harmonic_number(x,2)
+            harmonic_number(x, 2)
+        """
+        return BuiltinFunction.__call__(self, z, m, **kwds)
+
+    def _eval_(self, z, m):
+        """
+        EXAMPLES::
+
+            sage: harmonic_number(x,0)
+            x
+            sage: harmonic_number(x,1)
+            harmonic_number(x)
+            sage: harmonic_number(5)
+            137/60
+            sage: harmonic_number(3,3)
+            251/216
+            sage: harmonic_number(3,3).n() # this goes from rational to float
+            1.16203703703704
+            sage: harmonic_number(3,3.) # the following uses zeta functions
+            1.16203703703704
+            sage: harmonic_number(3.,3)
+            zeta(3) - 0.0400198661225573
+            sage: harmonic_number(0.1,5)
+            zeta(5) - 0.650300133161038
+            sage: harmonic_number(0.1,5).n()
+            0.386627621982332
+            sage: harmonic_number(3,5/2)
+            1/27*sqrt(3) + 1/8*sqrt(2) + 1
+        """
+        if m == 0:
+            return z
+        elif m == 1:
+            return harmonic_m1._eval_(z)
+
+        if z in ZZ and z >= 0:
+            return sum(1/(k**m) for k in xrange(1,z+1))
+
+    def _evalf_(self, z, m, parent=None, algorithm=None):
+        """
+        EXAMPLES::
+
+            sage: harmonic_number(3.,3)
+            zeta(3) - 0.0400198661225573
+            sage: harmonic_number(3.,3.)
+            1.16203703703704
+            sage: harmonic_number(3,3).n(200)
+            1.16203703703703703703703...
+            sage: harmonic_number(5,I).n()
+            2.36889632899995 - 3.51181956521611*I
+        """
+        if m == 0:
+            if parent is None:
+                return z
+            return parent(z)
+        elif m == 1:
+            return harmonic_m1._evalf_(z, parent, algorithm)
+
+        from sage.functions.transcendental import zeta, hurwitz_zeta
+        return zeta(m) - hurwitz_zeta(m,z+1)
+
+    def _maxima_init_evaled_(self, n, z):
+        """
+        EXAMPLES:
+
+            sage: maxima_calculus(harmonic_number(x,2))
+            gen_harmonic_number(2,_SAGE_VAR_x)
+            sage: maxima_calculus(harmonic_number(3,harmonic_number(x,3),hold=True))
+            1/3^gen_harmonic_number(3,_SAGE_VAR_x)+1/2^gen_harmonic_number(3,_SAGE_VAR_x)+1
+        """
+        if isinstance(n,str):
+            maxima_n=n
+        elif hasattr(n,'_maxima_init_'):
+            maxima_n=n._maxima_init_()
+        else:
+            maxima_n=str(n)
+        if isinstance(z,str):
+            maxima_z=z
+        elif hasattr(z,'_maxima_init_'):
+            maxima_z=z._maxima_init_()
+        else:
+            maxima_z=str(z)
+        return "gen_harmonic_number(%s,%s)" % (maxima_z, maxima_n) # swap arguments
+
+    def _derivative_(self, n, m, diff_param=None):
+        """
+        The derivative of `H_{n,m}`.
+
+        EXAMPLES::
+
+            sage: k,m,n = var('k,m,n')
+            sage: sum(1/k, k, 1, x).diff(x)
+            1/6*pi^2 - harmonic_number(x, 2)
+            sage: harmonic_number(x, 1).diff(x)
+            1/6*pi^2 - harmonic_number(x, 2)
+            sage: harmonic_number(n, m).diff(n)
+            -m*(harmonic_number(n, m + 1) - zeta(m + 1))
+            sage: harmonic_number(n, m).diff(m)
+            Traceback (most recent call last):
+            ...
+            ValueError: cannot differentiate harmonic_number in the second parameter
+        """
+        from sage.functions.transcendental import zeta
+        if diff_param == 1:
+            raise ValueError("cannot differentiate harmonic_number in the second parameter")
+        if m==1:
+            return harmonic_m1(n).diff()
+        else:
+            return m*(zeta(m+1) - harmonic_number(n, m+1))
+
+    def _print_(self, z, m):
+        """
+        EXAMPLES::
+
+            sage: harmonic_number(x)
+            harmonic_number(x)
+            sage: harmonic_number(x,2)
+            harmonic_number(x, 2)
+        """
+        if m == 1:
+            return "harmonic_number(%s)" % z
+        else:
+            return "harmonic_number(%s, %s)" % (z, m)
+
+    def _print_latex_(self, z, m):
+        """
+        EXAMPLES::
+
+            sage: latex(harmonic_number(x))
+            H_{x}
+            sage: latex(harmonic_number(x,2))
+            H_{{x},{2}}
+        """
+        if m == 1:
+            return r"H_{%s}" % z
+        else:
+            return r"H_{{%s},{%s}}" % (z, m)
+
+harmonic_number = Function_harmonic_number_generalized()
+
+def _swap_harmonic(a,b): return harmonic_number(b,a)
+
+from sage.symbolic.pynac import register_symbol
+
+register_symbol(_swap_harmonic,{'maxima':'gen_harmonic_number'})
+register_symbol(_swap_harmonic,{'maple':'harmonic'})
+
+class Function_harmonic_number(BuiltinFunction):
+    r"""
+    Harmonic number function, defined by:
+
+    .. math::
+
+        H_{n}=H_{n,1}=\sum_{k=1}^n\frac1k
+
+        H_{s}=\int_0^1\frac{1-x^s}{1-x}
+
+    See the docstring for :meth:`Function_harmonic_number_generalized`.
+
+    This class exists as callback for ``harmonic_number`` returned by Maxima.
+    """
+
+    def __init__(self):
+        r"""
+        EXAMPLES::
+
+            sage: k=var('k')
+            sage: loads(dumps(sum(1/k,k,1,x)))
+            harmonic_number(x)
+        """
+        BuiltinFunction.__init__(self, "harmonic_number", nargs=1,
+                                 conversions={'mathematica':'HarmonicNumber',
+                                              'maple':'harmonic',
+                                              'maxima':'harmonic_number'})
+
+    def _eval_(self, z, **kwds):
+        """
+        EXAMPLES::
+
+            sage: harmonic_number(0)
+            0
+            sage: harmonic_number(1)
+            1
+            sage: harmonic_number(20)
+            55835135/15519504
+            sage: harmonic_number(5/2)
+            -2*log(2) + 46/15
+            sage: harmonic_number(2*x)
+            harmonic_number(2*x)
+        """
+        if z in ZZ:
+            if z == 0:
+                return Integer(0)
+            elif z == 1:
+                return Integer(1)
+            elif z > 1:
+                import sage.libs.flint.arith as flint_arith
+                return flint_arith.harmonic_number(z)
+        elif z in QQ:
+            from sage.functions.other import psi1
+            return psi1(z+1) - psi1(1)
+
+    def _evalf_(self, z, parent=None, algorithm='mpmath'):
+        """
+        EXAMPLES::
+
+            sage: harmonic_number(20).n() # this goes from rational to float
+            3.59773965714368
+            sage: harmonic_number(20).n(200)
+            3.59773965714368191148376906...
+            sage: harmonic_number(20.) # this computes the integral with mpmath
+            3.59773965714368
+            sage: harmonic_number(1.0*I)
+            0.671865985524010 + 1.07667404746858*I
+        """
+        from sage.libs.mpmath import utils as mpmath_utils
+        import mpmath
+        return mpmath_utils.call(mpmath.harmonic, z, parent=parent)
+
+    def _derivative_(self, z, diff_param=None):
+        """
+        The derivative of `H_x`.
+
+        EXAMPLES::
+
+            sage: k=var('k')
+            sage: sum(1/k,k,1,x).diff(x)
+            1/6*pi^2 - harmonic_number(x, 2)
+        """
+        from sage.functions.transcendental import zeta
+        return zeta(2)-harmonic_number(z,2)
+
+    def _print_latex_(self, z):
+        """
+        EXAMPLES::
+
+            sage: k=var('k')
+            sage: latex(sum(1/k,k,1,x))
+            H_{x}
+        """
+        return r"H_{%s}" % z
+
+harmonic_m1 = Function_harmonic_number()
diff --git a/src/sage/interfaces/maxima_lib.py b/src/sage/interfaces/maxima_lib.py
index 1f488b99034..f3a6862ead9 100644
--- a/src/sage/interfaces/maxima_lib.py
+++ b/src/sage/interfaces/maxima_lib.py
@@ -1255,6 +1255,7 @@ def sage_rat(x,y):
 max_array = EclObject("ARRAY")
 mdiff = EclObject("%DERIVATIVE")
 max_lambert_w = sage_op_dict[sage.functions.log.lambert_w]
+max_harmo = EclObject("$GEN_HARMONIC_NUMBER")
 
 def mrat_to_sage(expr):
     r"""
@@ -1436,6 +1437,20 @@ def dummy_integrate(expr):
         return sage.symbolic.integration.integral.indefinite_integral(*args,
                                                                   hold=True)
 
+def max_harmonic_to_sage(expr):
+    """
+    EXAMPLES::
+
+        sage: from sage.interfaces.maxima_lib import maxima_lib, max_to_sr
+        sage: c=maxima_lib(harmonic_number(x,2))
+        sage: c.ecl()
+        <ECL: (($GEN_HARMONIC_NUMBER SIMP) 2 |$_SAGE_VAR_x|)>
+        sage: max_to_sr(c.ecl())
+        harmonic_number(x, 2)
+    """
+    return sage.functions.log.harmonic_number(max_to_sr(caddr(expr)),
+                                              max_to_sr(cadr(expr)))
+
 ## The dictionaries
 special_max_to_sage={
     mrat : mrat_to_sage,
@@ -1443,7 +1458,8 @@ def dummy_integrate(expr):
     mdiff : mdiff_to_sage,
     EclObject("%INTEGRATE") : dummy_integrate,
     max_at : max_at_to_sage,
-    mlist : mlist_to_sage
+    mlist : mlist_to_sage,
+    max_harmo : max_harmonic_to_sage
 }
 
 special_sage_to_max={
@@ -1451,6 +1467,7 @@ def dummy_integrate(expr):
     sage.functions.other.psi1 : lambda X : [[mqapply],[[max_psi, max_array],0],X],
     sage.functions.other.psi2 : lambda N,X : [[mqapply],[[max_psi, max_array],N],X],
     sage.functions.log.lambert_w : lambda N,X : [[max_lambert_w], X] if N==EclObject(0) else [[mqapply],[[max_lambert_w, max_array],N],X],
+    sage.functions.log.harmonic_number : lambda N,X : [[max_harmo],X,N],
     sage.functions.hypergeometric.hypergeometric : lambda A, B, X : [[mqapply],[[max_hyper, max_array],lisp_length(A.cdr()),lisp_length(B.cdr())],A,B,X]
 }
 
diff --git a/src/sage/libs/flint/arith.pxd b/src/sage/libs/flint/arith.pxd
index 6220d96ae59..4cd331576a6 100644
--- a/src/sage/libs/flint/arith.pxd
+++ b/src/sage/libs/flint/arith.pxd
@@ -8,3 +8,4 @@ cdef extern from "flint/arith.h":
     void arith_euler_number ( fmpz_t res , ulong n )
     void arith_number_of_partitions(fmpz_t x, ulong n)
     void arith_dedekind_sum(fmpq_t, fmpz_t, fmpz_t)
+    void arith_harmonic_number(fmpq_t, unsigned long n)
diff --git a/src/sage/libs/flint/arith.pyx b/src/sage/libs/flint/arith.pyx
index 934a8f8a2c1..f5156008232 100644
--- a/src/sage/libs/flint/arith.pyx
+++ b/src/sage/libs/flint/arith.pyx
@@ -17,6 +17,7 @@ include "cysignals/signals.pxi"
 from .fmpz cimport *
 from .fmpq cimport *
 
+
 from sage.rings.integer cimport Integer
 from sage.rings.rational cimport Rational
 
@@ -207,3 +208,29 @@ def dedekind_sum(p, q):
 
     return s
 
+def harmonic_number(unsigned long n):
+    """
+    Returns the harmonic number ``H_n``.
+
+    EXAMPLES::
+
+        sage: from sage.libs.flint.arith import harmonic_number
+        sage: n = 500 + randint(0,500)
+        sage: bool( sum(1/k for k in range(1,n+1)) == harmonic_number(n) )
+        True
+    """
+    s = Rational(0)
+    cdef fmpq_t s_fmpq
+
+    fmpq_init(s_fmpq)
+
+    sig_on()
+    arith_harmonic_number(s_fmpq, n)
+
+    fmpq_get_mpq((<Rational>s).value, s_fmpq)
+    sig_off()
+
+    fmpq_clear(s_fmpq)
+
+    return s
+