def test_issue_14177():
    n = Symbol('n', positive=True, integer=True)

    assert zeta(2*n) == (-1)**(n + 1)*2**(2*n - 1)*pi**(2*n)*bernoulli(2*n)/factorial(2*n)
    assert zeta(-n) == (-1)**(-n)*bernoulli(n + 1)/(n + 1)

    n = Symbol('n')

    assert zeta(2*n) == zeta(2*n) # As sign of z (= 2*n) is not determined

    def _eval_aseries(self, n, args0, x, logx):
        from sympy import Order
        if args0[1] != oo or not \
                (self.args[0].is_Integer and self.args[0].is_nonnegative):
            return super(polygamma, self)._eval_aseries(n, args0, x, logx)
        z = self.args[1]
        N = self.args[0]

        if N == 0:
            # digamma function series
            # Abramowitz & Stegun, p. 259, 6.3.18
            r = log(z) - 1/(2*z)
            o = None
            if n < 2:
                o = Order(1/z, x)
            else:
                m = ceiling((n + 1)//2)
                l = [bernoulli(2*k) / (2*k*z**(2*k)) for k in range(1, m)]
                r -= Add(*l)
                o = Order(1/z**(2*m), x)
            return r._eval_nseries(x, n, logx) + o
        else:
            # proper polygamma function
            # Abramowitz & Stegun, p. 260, 6.4.10
            # We return terms to order higher than O(x**n) on purpose
            # -- otherwise we would not be able to return any terms for
            #    quite a long time!
            fac = gamma(N)
            e0 = fac + N*fac/(2*z)
            m = ceiling((n + 1)//2)
            for k in range(1, m):
                fac = fac*(2*k + N - 1)*(2*k + N - 2) / ((2*k)*(2*k - 1))
                e0 += bernoulli(2*k)*fac/z**(2*k)
            o = Order(1/z**(2*m), x)
            if n == 0:
                o = Order(1/z, x)
            elif n == 1:
                o = Order(1/z**2, x)
            r = e0._eval_nseries(z, n, logx) + o
            return (-1 * (-1/z)**N * r)._eval_nseries(x, n, logx)

    def eval(cls, z, a_=None):
        if a_ is None:
            z, a = list(map(sympify, (z, 1)))
        else:
            z, a = list(map(sympify, (z, a_)))

        if a.is_Number:
            if a is S.NaN:
                return S.NaN
            elif a is S.One and a_ is not None:
                return cls(z)
            # TODO Should a == 0 return S.NaN as well?

        if z.is_Number:
            if z is S.NaN:
                return S.NaN
            elif z is S.Infinity:
                return S.One
            elif z is S.Zero:
                if a.is_negative:
                    return S.Half - a - 1
                else:
                    return S.Half - a
            elif z is S.One:
                return S.ComplexInfinity
            elif z.is_Integer:
                if a.is_Integer:
                    if z.is_negative:
                        zeta = (-1)**z * bernoulli(-z + 1)/(-z + 1)
                    elif z.is_even:
                        B, F = bernoulli(z), factorial(z)
                        zeta = 2**(z - 1) * abs(B) * pi**z / F
                    else:
                        return

                    if a.is_negative:
                        return zeta + harmonic(abs(a), z)
                    else:
                        return zeta - harmonic(a - 1, z)

 def _eval_aseries(self, n, args0, x, logx):
     if args0[0] != oo:
         return super(loggamma, self)._eval_aseries(n, args0, x, logx)
     z = self.args[0]
     m = min(n, C.ceiling((n + S(1)) / 2))
     r = log(z) * (z - S(1) / 2) - z + log(2 * pi) / 2
     l = [bernoulli(2 * k) / (2 * k * (2 * k - 1) * z ** (2 * k - 1)) for k in range(1, m)]
     o = None
     if m == 0:
         o = C.Order(1, x)
     else:
         o = C.Order(1 / z ** (2 * m - 1), x)
     # It is very inefficient to first add the order and then do the nseries
     return (r + Add(*l))._eval_nseries(x, n, logx) + o

def test_sympy__functions__combinatorial__numbers__bernoulli():
    from sympy.functions.combinatorial.numbers import bernoulli
    assert _test_args(bernoulli(x))