numbers.py

"""
This module implements some special functions that commonly appear in
combinatorial contexts (e.g. in power series); in particular,
sequences of rational numbers such as Bernoulli and Fibonacci numbers.

Factorials, binomial coefficients and related functions are located in
the separate 'factorials' module.
"""

from __future__ import annotations

import collections

from mpmath import bernfrac, mp, workprec
from mpmath.libmp import ifib as _ifib

from ...core import (Add, Dummy, E, Expr, Function, GoldenRatio, Integer,
                     Rational, cacheit, expand_mul, nan, oo, pi)
from ...core.compatibility import as_int
from ...utilities.memoization import recurrence_memo
from ..elementary.exponential import log
from ..elementary.integers import floor
from ..elementary.trigonometric import cos, cot, sin
from .factorials import binomial, factorial


def _product(a, b):
    p = 1
    for k in range(a, b + 1):
        p *= k
    return p


# Dummy symbol used for computing polynomial sequences
_sym = Dummy('_for_recurrence_memo')


############################################################################
#                                                                          #
#                           Fibonacci numbers                              #
#                                                                          #
############################################################################

class fibonacci(Function):
    r"""
    Fibonacci numbers / Fibonacci polynomials

    The Fibonacci numbers are the integer sequence defined by the
    initial terms F_0 = 0, F_1 = 1 and the two-term recurrence
    relation F_n = F_{n-1} + F_{n-2}.  This definition
    extended to arbitrary real and complex arguments using
    the formula

    .. math :: F_z = \frac{\phi^z - \cos(\pi z) \phi^{-z}}{\sqrt 5}

    The Fibonacci polynomials are defined by F_1(x) = 1,
    F_2(x) = x, and F_n(x) = x*F_{n-1}(x) + F_{n-2}(x) for n > 2.
    For all positive integers n, F_n(1) = F_n.

    * fibonacci(n) gives the nth Fibonacci number, F_n
    * fibonacci(n, x) gives the nth Fibonacci polynomial in x, F_n(x)

    Examples
    ========

    >>> [fibonacci(x) for x in range(11)]
    [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
    >>> fibonacci(5, Symbol('t'))
    t**4 + 3*t**2 + 1

    References
    ==========

    * https://en.wikipedia.org/wiki/Fibonacci_number
    * https://mathworld.wolfram.com/FibonacciNumber.html

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.catalan
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.harmonic
    diofant.functions.combinatorial.numbers.lucas

    """

    @staticmethod
    def _fib(n):
        return _ifib(n)

    @staticmethod
    @recurrence_memo([None, Integer(1), _sym])
    def _fibpoly(n, prev):
        return (prev[-2] + _sym*prev[-1]).expand()

    @classmethod
    def eval(cls, n, sym=None):
        if n.is_Integer:
            n = int(n)
            if sym is None:
                if n < 0:
                    return (-1)**(n + 1) * fibonacci(-n)
                return Integer(cls._fib(n))
            if n < 1:
                raise ValueError('Fibonacci polynomials are defined '
                                 'only for positive integer indices.')
            return cls._fibpoly(n).subs({_sym: sym})

    def _eval_rewrite_as_sqrt(self, n, sym=None, **kwargs):
        from .. import sqrt
        if sym is None:
            return (GoldenRatio**n - cos(pi*n)/GoldenRatio**n)/sqrt(5)
    _eval_rewrite_as_tractable = _eval_rewrite_as_sqrt


class lucas(Function):
    """
    Lucas numbers

    Lucas numbers satisfy a recurrence relation similar to that of
    the Fibonacci sequence, in which each term is the sum of the
    preceding two. They are generated by choosing the initial
    values L_0 = 2 and L_1 = 1.

    * lucas(n) gives the nth Lucas number

    Examples
    ========

    >>> [lucas(x) for x in range(11)]
    [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123]

    References
    ==========

    * https://en.wikipedia.org/wiki/Lucas_number
    * https://mathworld.wolfram.com/LucasNumber.html

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.catalan
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.harmonic

    """

    @classmethod
    def eval(cls, n):
        if n.is_Integer:
            return fibonacci(n + 1) + fibonacci(n - 1)


############################################################################
#                                                                          #
#                           Bernoulli numbers                              #
#                                                                          #
############################################################################

class bernoulli(Function):
    r"""
    Bernoulli numbers / Bernoulli polynomials

    The Bernoulli numbers are a sequence of rational numbers
    defined by B_0 = 1 and the recursive relation (n > 0)::

                n
               ___
              \      / n + 1 \
          0 =  )     |       | * B .
              /___   \   k   /    k
              k = 0

    They are also commonly defined by their exponential generating
    function, which is x/(exp(x) - 1). For odd indices > 1, the
    Bernoulli numbers are zero.

    The Bernoulli polynomials satisfy the analogous formula::

                    n
                   ___
                  \      / n \         n-k
          B (x) =  )     |   | * B  * x   .
           n      /___   \ k /    k
                  k = 0

    Bernoulli numbers and Bernoulli polynomials are related as
    B_n(0) = B_n.

    We compute Bernoulli numbers using Ramanujan's formula::

                                   / n + 3 \
          B   =  (A(n) - S(n))  /  |       |
           n                       \   n   /

    where A(n) = (n+3)/3 when n = 0 or 2 (mod 6), A(n) = -(n+3)/6
    when n = 4 (mod 6), and::

                 [n/6]
                  ___
                 \      /  n + 3  \
          S(n) =  )     |         | * B
                 /___   \ n - 6*k /    n-6*k
                 k = 1

    This formula is similar to the sum given in the definition, but
    cuts 2/3 of the terms. For Bernoulli polynomials, we use the
    formula in the definition.

    * bernoulli(n) gives the nth Bernoulli number, B_n
    * bernoulli(n, x) gives the nth Bernoulli polynomial in x, B_n(x)

    Examples
    ========

    >>> [bernoulli(n) for n in range(11)]
    [1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
    >>> bernoulli(1000001)
    0

    References
    ==========

    * https://en.wikipedia.org/wiki/Bernoulli_number
    * https://en.wikipedia.org/wiki/Bernoulli_polynomial
    * https://mathworld.wolfram.com/BernoulliNumber.html
    * https://mathworld.wolfram.com/BernoulliPolynomial.html

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.catalan
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.harmonic
    diofant.functions.combinatorial.numbers.lucas

    """

    # Calculates B_n for positive even n
    @staticmethod
    def _calc_bernoulli(n):
        s = 0
        a = int(binomial(n + 3, n - 6))
        for j in range(1, n//6 + 1):
            s += a * bernoulli(n - 6*j)
            # Avoid computing each binomial coefficient from scratch
            a *= _product(n - 6 - 6*j + 1, n - 6*j)
            a //= _product(6*j + 4, 6*j + 9)
        if n % 6 == 4:
            s = -Rational(n + 3, 6) - s
        else:
            s = Rational(n + 3, 3) - s
        return s / binomial(n + 3, n)

    # We implement a specialized memoization scheme to handle each
    # case modulo 6 separately
    _cache = {0: Integer(1), 2: Rational(1, 6), 4: Rational(-1, 30)}
    _highest = {0: 0, 2: 2, 4: 4}

    @classmethod
    def eval(cls, n, sym=None):
        if n.is_Number:
            if n.is_Integer and n.is_nonnegative:
                if n == 0:
                    return Integer(1)
                if n == 1:
                    if sym is None:
                        return -Rational(1, 2)
                    return sym - Rational(1, 2)
                # Bernoulli numbers
                if sym is None:
                    if n.is_odd:
                        return Integer(0)
                    n = int(n)
                    # Use mpmath for enormous Bernoulli numbers
                    if n > 500:
                        p, q = bernfrac(n)
                        return Rational(int(p), int(q))
                    case = n % 6
                    highest_cached = cls._highest[case]
                    if n <= highest_cached:
                        return cls._cache[n]
                    # To avoid excessive recursion when, say, bernoulli(1000) is
                    # requested, calculate and cache the entire sequence ... B_988,
                    # B_994, B_1000 in increasing order
                    for i in range(highest_cached + 6, n + 6, 6):
                        b = cls._calc_bernoulli(i)
                        cls._cache[i] = b
                        cls._highest[case] = i
                    return b
                # Bernoulli polynomials
                n, result = int(n), []
                for k in range(n + 1):
                    result.append(binomial(n, k)*cls(k)*sym**(n - k))
                return Add(*result)
            raise ValueError('Bernoulli numbers are defined only'
                             ' for nonnegative integer indices.')

        if sym is None:
            if n.is_odd and (n - 1).is_positive:
                return Integer(0)


############################################################################
#                                                                          #
#                             Bell numbers                                 #
#                                                                          #
############################################################################

class bell(Function):
    r"""
    Bell numbers / Bell polynomials

    The Bell numbers satisfy `B_0 = 1` and

    .. math:: B_n = \sum_{k=0}^{n-1} \binom{n-1}{k} B_k.

    They are also given by:

    .. math:: B_n = \frac{1}{e} \sum_{k=0}^{\infty} \frac{k^n}{k!}.

    The Bell polynomials are given by `B_0(x) = 1` and

    .. math:: B_n(x) = x \sum_{k=1}^{n-1} \binom{n-1}{k-1} B_{k-1}(x).

    The second kind of Bell polynomials (are sometimes called "partial" Bell
    polynomials or incomplete Bell polynomials) are defined as

    .. math:: B_{n,k}(x_1, x_2,\dotsc x_{n-k+1}) =
            \sum_{j_1+j_2+j_2+\dotsb=k \atop j_1+2j_2+3j_2+\dotsb=n}
                \frac{n!}{j_1!j_2!\dotsb j_{n-k+1}!}
                \left(\frac{x_1}{1!} \right)^{j_1}
                \left(\frac{x_2}{2!} \right)^{j_2} \dotsb
                \left(\frac{x_{n-k+1}}{(n-k+1)!} \right) ^{j_{n-k+1}}.

    * bell(n) gives the `n^{th}` Bell number, `B_n`.
    * bell(n, x) gives the `n^{th}` Bell polynomial, `B_n(x)`.
    * bell(n, k, (x1, x2, ...)) gives Bell polynomials of the second kind,
      `B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1})`.

    Notes
    =====

    Not to be confused with Bernoulli numbers and Bernoulli polynomials,
    which use the same notation.

    Examples
    ========

    >>> [bell(n) for n in range(11)]
    [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975]
    >>> bell(30)
    846749014511809332450147
    >>> bell(4, Symbol('t'))
    t**4 + 6*t**3 + 7*t**2 + t
    >>> bell(6, 2, symbols('x:6')[1:])
    6*x1*x5 + 15*x2*x4 + 10*x3**2

    References
    ==========

    * https://en.wikipedia.org/wiki/Bell_number
    * https://mathworld.wolfram.com/BellNumber.html
    * https://mathworld.wolfram.com/BellPolynomial.html

    See Also
    ========

    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.catalan
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.harmonic
    diofant.functions.combinatorial.numbers.lucas

    """

    @staticmethod
    @recurrence_memo([1, 1])
    def _bell(n, prev):
        s = 1
        a = 1
        for k in range(1, n):
            a = a * (n - k) // k
            s += a * prev[k]
        return s

    @staticmethod
    @recurrence_memo([Integer(1), _sym])
    def _bell_poly(n, prev):
        s = 1
        a = 1
        for k in range(2, n + 1):
            a = a * (n - k + 1) // (k - 1)
            s += a * prev[k - 1]
        return expand_mul(_sym * s)

    @staticmethod
    def _bell_incomplete_poly(n, k, symbols):
        r"""
        The second kind of Bell polynomials (incomplete Bell polynomials).

        Calculated by recurrence formula:

        .. math:: B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1}) =
                \sum_{m=1}^{n-k+1}
                \x_m \binom{n-1}{m-1} B_{n-m,k-1}(x_1, x_2, \dotsc, x_{n-m-k})

        where
            B_{0,0} = 1;
            B_{n,0} = 0; for n>=1
            B_{0,k} = 0; for k>=1

        """
        if (n == 0) and (k == 0):
            return Integer(1)
        if (n == 0) or (k == 0):
            return Integer(0)
        s = Integer(0)
        a = Integer(1)
        for m in range(1, n - k + 2):
            s += a * bell._bell_incomplete_poly(
                n - m, k - 1, symbols) * symbols[m - 1]
            a = a * (n - m) / m
        return expand_mul(s)

    @classmethod
    def eval(cls, n, k_sym=None, symbols=None):
        if n.is_Integer and n.is_nonnegative:
            if k_sym is None:
                return Integer(cls._bell(int(n)))
            if symbols is None:
                return cls._bell_poly(int(n)).subs({_sym: k_sym})
            r = cls._bell_incomplete_poly(int(n), int(k_sym), symbols)
            return r

    def _eval_rewrite_as_Sum(self, n, k_sym=None, symbols=None):
        from ...concrete import Sum
        if (k_sym is not None) or (symbols is not None):
            return self

        # Dobinski's formula
        if not n.is_nonnegative:
            return self
        k = Dummy('k', integer=True, nonnegative=True)
        return 1 / E * Sum(k**n / factorial(k), (k, 0, oo))

############################################################################
#                                                                          #
#                           Harmonic numbers                               #
#                                                                          #
############################################################################


class harmonic(Function):
    r"""
    Harmonic numbers

    The nth harmonic number is given by `\operatorname{H}_{n} =
    1 + \frac{1}{2} + \frac{1}{3} + \ldots + \frac{1}{n}`.

    More generally:

    .. math:: \operatorname{H}_{n,m} = \sum_{k=1}^{n} \frac{1}{k^m}

    As `n \rightarrow \infty`, `\operatorname{H}_{n,m} \rightarrow \zeta(m)`,
    the Riemann zeta function.

    * ``harmonic(n)`` gives the nth harmonic number, `\operatorname{H}_n`

    * ``harmonic(n, m)`` gives the nth generalized harmonic number
      of order `m`, `\operatorname{H}_{n,m}`, where
      ``harmonic(n) == harmonic(n, 1)``

    Examples
    ========

    >>> [harmonic(n) for n in range(6)]
    [0, 1, 3/2, 11/6, 25/12, 137/60]
    >>> [harmonic(n, 2) for n in range(6)]
    [0, 1, 5/4, 49/36, 205/144, 5269/3600]
    >>> harmonic(oo, 2)
    pi**2/6

    >>> harmonic(n).rewrite(Sum)
    Sum(1/_k, (_k, 1, n))

    We can evaluate harmonic numbers for all integral and positive
    rational arguments:

    >>> harmonic(8)
    761/280
    >>> harmonic(11)
    83711/27720

    >>> H = harmonic(Rational(1, 3))
    >>> H
    harmonic(1/3)
    >>> He = expand_func(H)
    >>> He
    -log(6) - sqrt(3)*pi/6 + 2*Sum(log(sin(pi*_k/3))*cos(2*pi*_k/3), (_k, 1, 1))
                           + 3*Sum(1/(3*_k + 1), (_k, 0, 0))
    >>> He.doit()
    -log(6) - sqrt(3)*pi/6 - log(sqrt(3)/2) + 3
    >>> H = harmonic(Rational(25, 7))
    >>> He = simplify(expand_func(H).doit())
    >>> He
    log(sin(pi/7)**(-2*cos(pi/7))*sin(2*pi/7)**(2*cos(16*pi/7))*cos(pi/14)**(-2*sin(pi/14))/14)
    + pi*tan(pi/14)/2 + 30247/9900
    >>> He.evalf(40)
    1.983697455232980674869851942390639915940
    >>> harmonic(Rational(25, 7)).evalf(40)
    1.983697455232980674869851942390639915940

    We can rewrite harmonic numbers in terms of polygamma functions:

    >>> harmonic(n).rewrite(digamma)
    polygamma(0, n + 1) + EulerGamma

    >>> harmonic(n).rewrite(polygamma)
    polygamma(0, n + 1) + EulerGamma

    >>> harmonic(n, 3).rewrite(polygamma)
    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)

    Integer offsets in the argument can be pulled out:

    >>> expand_func(harmonic(n+4))
    harmonic(n) + 1/(n + 4) + 1/(n + 3) + 1/(n + 2) + 1/(n + 1)

    >>> expand_func(harmonic(n-4))
    harmonic(n) - 1/(n - 1) - 1/(n - 2) - 1/(n - 3) - 1/n

    Some limits can be computed as well:

    >>> limit(harmonic(n), n, oo)
    oo

    >>> limit(harmonic(n, 2), n, oo)
    pi**2/6

    >>> limit(harmonic(n, 3), n, oo)
    -polygamma(2, 1)/2

    However we can not compute the general relation yet:

    >>> limit(harmonic(n, m), n, oo)
    harmonic(oo, m)

    which equals ``zeta(m)`` for ``m > 1``.

    References
    ==========

    * https://en.wikipedia.org/wiki/Harmonic_number
    * http://functions.wolfram.com/GammaBetaErf/HarmonicNumber/
    * http://functions.wolfram.com/GammaBetaErf/HarmonicNumber2/

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.catalan
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.lucas

    """

    # Generate one memoized Harmonic number-generating function for each
    # order and store it in a dictionary
    _functions: dict[Integer, collections.abc.Callable[[int], Rational]] = {}

    @classmethod
    def eval(cls, n, m=None):
        from .. import zeta
        if m == 1:
            return cls(n)
        if m is None:
            m = Integer(1)

        if m.is_zero:
            return n

        if n is oo:
            if m.is_negative:
                return nan
            if (m - 1).is_nonpositive and m.is_nonnegative:
                return oo
            if (m - 1).is_positive:
                return zeta(m)

        if n.is_Integer and n.is_nonnegative and m.is_Integer:
            if n == 0:
                return Integer(0)
            if m not in cls._functions:
                @recurrence_memo([0])
                def f(n, prev):
                    return prev[-1] + 1/n**m
                cls._functions[m] = f
            return cls._functions[m](int(n))

    def _eval_rewrite_as_polygamma(self, n, m=1):
        from .. import polygamma
        return (-1)**m/factorial(m - 1) * (polygamma(m - 1, 1) - polygamma(m - 1, n + 1))

    def _eval_rewrite_as_digamma(self, n, m=1):
        from .. import polygamma
        return self.rewrite(polygamma)

    def _eval_rewrite_as_trigamma(self, n, m=1):
        from .. import polygamma
        return self.rewrite(polygamma)

    def _eval_rewrite_as_Sum(self, n, m=None):
        from ...concrete import Sum
        k = Dummy('k')
        if m is None:
            m = Integer(1)
        return Sum(k**(-m), (k, 1, n))

    def _eval_expand_func(self, **hints):
        from ...concrete import Sum
        n = self.args[0]
        m = self.args[1] if len(self.args) == 2 else 1

        if m == 1:
            if n.is_Add:
                off = n.args[0]
                nnew = n - off
                if off.is_Integer and off.is_positive:
                    result = [1/(nnew + i) for i in range(off, 0, -1)] + [harmonic(nnew)]
                    return Add(*result)
                if off.is_Integer and off.is_negative:
                    result = [-1/(nnew + i) for i in range(0, off, -1)] + [harmonic(nnew)]
                    return Add(*result)

            if n.is_Rational:
                # Expansions for harmonic numbers at general rational arguments (u + p/q)
                # Split n as u + p/q with p < q
                p, q = n.as_numer_denom()
                u = p // q
                p = p - u * q
                if u.is_nonnegative and p.is_positive and q.is_positive and p < q:
                    k = Dummy('k')
                    t1 = q * Sum(1 / (q * k + p), (k, 0, u))
                    t2 = 2 * Sum(cos((2 * pi * p * k) / q) *
                                 log(sin((pi * k) / q)),
                                 (k, 1, floor(Rational(q - 1, 2))))
                    t3 = (pi / 2) * cot((pi * p) / q) + log(2 * q)
                    return t1 + t2 - t3

        return self

    def _eval_rewrite_as_tractable(self, n, m=1, **kwargs):
        from .. import polygamma
        return self.rewrite(polygamma).rewrite('tractable')

    def _eval_evalf(self, prec):
        from .. import polygamma
        if all(i.is_number for i in self.args):
            return self.rewrite(polygamma)._eval_evalf(prec)


############################################################################
#                                                                          #
#                           Euler numbers                                  #
#                                                                          #
############################################################################


class euler(Function):
    r"""
    Euler numbers

    The euler numbers are given by::

                  2*n+1   k
                   ___   ___            j          2*n+1
                  \     \     / k \ (-1)  * (k-2*j)
          E   = I  )     )    |   | --------------------
           2n     /___  /___  \ j /      k    k
                  k = 1 j = 0           2  * I  * k

          E     = 0
           2n+1

    * euler(n) gives the n-th Euler number, E_n

    Examples
    ========

    >>> from diofant.functions import euler

    >>> [euler(n) for n in range(10)]
    [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
    >>> euler(n+2*n)
    euler(3*n)

    References
    ==========

    * https://en.wikipedia.org/wiki/Euler_numbers
    * https://mathworld.wolfram.com/EulerNumber.html
    * https://en.wikipedia.org/wiki/Alternating_permutation
    * https://mathworld.wolfram.com/AlternatingPermutation.html

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.harmonic
    diofant.functions.combinatorial.numbers.lucas

    """

    @classmethod
    def eval(cls, m):
        if m.is_odd:
            return Integer(0)
        if m.is_Integer and m.is_nonnegative:
            m = m._to_mpmath(mp.prec)
            res = mp.eulernum(m, exact=True)
            return Integer(res)

    def _eval_evalf(self, prec):
        m = self.args[0]

        if m.is_Integer and m.is_nonnegative:
            m = m._to_mpmath(prec)
            with workprec(prec):
                res = mp.eulernum(m)
            return Expr._from_mpmath(res, prec)

############################################################################
#                                                                          #
#                           Catalan numbers                                #
#                                                                          #
############################################################################


class catalan(Function):
    r"""
    Catalan numbers

    The n-th catalan number is given by::

                 1   / 2*n \
          C  = ----- |     |
           n   n + 1 \  n  /

    * catalan(n) gives the n-th Catalan number, C_n

    Examples
    ========

    >>> [catalan(i) for i in range(1, 10)]
    [1, 2, 5, 14, 42, 132, 429, 1430, 4862]

    >>> catalan(n)
    catalan(n)

    Catalan numbers can be transformed into several other, identical
    expressions involving other mathematical functions

    >>> catalan(n).rewrite(binomial)
    binomial(2*n, n)/(n + 1)

    >>> catalan(n).rewrite(gamma)
    4**n*gamma(n + 1/2)/(sqrt(pi)*gamma(n + 2))

    >>> catalan(n).rewrite(hyper)
    hyper((-n + 1, -n), (2,), 1)

    For some non-integer values of n we can get closed form
    expressions by rewriting in terms of gamma functions:

    >>> catalan(Rational(1, 2)).rewrite(gamma)
    8/(3*pi)

    We can differentiate the Catalan numbers C(n) interpreted as a
    continuous real function in n:

    >>> diff(catalan(n), n)
    (polygamma(0, n + 1/2) - polygamma(0, n + 2) + log(4))*catalan(n)

    As a more advanced example consider the following ratio
    between consecutive numbers:

    >>> combsimp((catalan(n + 1)/catalan(n)).rewrite(binomial))
    2*(2*n + 1)/(n + 2)

    The Catalan numbers can be generalized to complex numbers:

    >>> catalan(I).rewrite(gamma)
    4**I*gamma(1/2 + I)/(sqrt(pi)*gamma(2 + I))

    and evaluated with arbitrary precision:

    >>> catalan(I).evalf(20)
    0.39764993382373624267 - 0.020884341620842555705*I

    References
    ==========

    * https://en.wikipedia.org/wiki/Catalan_number
    * https://mathworld.wolfram.com/CatalanNumber.html
    * http://functions.wolfram.com/GammaBetaErf/CatalanNumber/
    * http://geometer.org/mathcircles/catalan.pdf

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.harmonic
    diofant.functions.combinatorial.numbers.lucas
    diofant.functions.combinatorial.factorials.binomial

    """

    @classmethod
    def eval(cls, n):
        from .. import gamma
        if (n.is_Integer and n.is_nonnegative) or \
           (n.is_noninteger and n.is_negative):
            return 4**n*gamma(n + Rational(1, 2))/(gamma(Rational(1, 2))*gamma(n + 2))

        if n.is_integer and n.is_negative:
            if (n + 1).is_negative:
                return Integer(0)
            return -Rational(1, 2)

    def fdiff(self, argindex=1):
        from .. import log, polygamma
        n = self.args[0]
        return catalan(n)*(polygamma(0, n + Rational(1, 2)) - polygamma(0, n + 2) + log(4))

    def _eval_rewrite_as_binomial(self, n):
        return binomial(2*n, n)/(n + 1)

    def _eval_rewrite_as_factorial(self, n):
        return factorial(2*n) / (factorial(n+1) * factorial(n))

    def _eval_rewrite_as_gamma(self, n):
        from .. import gamma

        # The gamma function allows to generalize Catalan numbers to complex n
        return 4**n*gamma(n + Rational(1, 2))/(gamma(Rational(1, 2))*gamma(n + 2))

    def _eval_rewrite_as_hyper(self, n):
        from .. import hyper
        return hyper([1 - n, -n], [2], 1)

    def _eval_rewrite_as_Product(self, n):
        from ...concrete import Product
        if not (n.is_integer and n.is_nonnegative):
            return self
        k = Dummy('k', integer=True, positive=True)
        return Product((n + k) / k, (k, 2, n))

    def _eval_evalf(self, prec):
        from .. import gamma
        if self.args[0].is_number:
            return self.rewrite(gamma)._eval_evalf(prec)


############################################################################
#                                                                          #
#                           Genocchi numbers                               #
#                                                                          #
############################################################################


class genocchi(Function):
    r"""
    Genocchi numbers

    The Genocchi numbers are a sequence of integers G_n that satisfy the
    relation::

                           oo
                         ____
                         \   `
                 2*t      \         n
                ------ =   \   G_n*t
                 t         /   ------
                e  + 1    /      n!
                         /___,
                         n = 1

    Examples
    ========

    >>> [genocchi(n) for n in range(1, 9)]
    [1, -1, 0, 1, 0, -3, 0, 17]
    >>> n = Symbol('n', integer=True, positive=True)
    >>> genocchi(2 * n + 1)
    0

    References
    ==========

    * https://en.wikipedia.org/wiki/Genocchi_number
    * https://mathworld.wolfram.com/GenocchiNumber.html

    See Also
    ========

    diofant.functions.combinatorial.numbers.bell
    diofant.functions.combinatorial.numbers.bernoulli
    diofant.functions.combinatorial.numbers.catalan
    diofant.functions.combinatorial.numbers.euler
    diofant.functions.combinatorial.numbers.fibonacci
    diofant.functions.combinatorial.numbers.harmonic
    diofant.functions.combinatorial.numbers.lucas

    """

    @classmethod
    def eval(cls, n):
        if n.is_Number:
            if (not n.is_Integer) or n.is_nonpositive:
                raise ValueError('Genocchi numbers are defined only for ' +
                                 'positive integers')
            return 2*(1 - 2**n)*bernoulli(n)

        if n.is_odd and (n - 1).is_positive:
            return Integer(0)

        if (n - 1).is_zero:
            return Integer(1)

    def _eval_rewrite_as_bernoulli(self, n):
        if n.is_integer and n.is_nonnegative:
            return 2*(1 - 2**n)*bernoulli(n)

    def _eval_is_negative(self):
        n = self.args[0]
        if n.is_integer and n.is_positive:
            if n.is_even:
                return (n/2).is_odd

    def _eval_is_odd(self):
        n = self.args[0]
        if n.is_integer and n.is_positive:
            if n.is_even:
                return True


@cacheit
def _stirling1(n, k):
    if n == k == 0:
        return Integer(1)
    if 0 in (n, k):
        return Integer(0)
    n1 = n - 1

    # some special values
    if n == k:
        return Integer(1)
    if k == 1:
        return factorial(n1)
    if k == n1:
        return binomial(n, 2)
    if k == n - 2:
        return (3*n - 1)*binomial(n, 3)/4
    if k == n - 3:
        return binomial(n, 2)*binomial(n, 4)

    # general recurrence
    return n1*_stirling1(n1, k) + _stirling1(n1, k - 1)


@cacheit
def _stirling2(n, k):
    if n == k == 0:
        return Integer(1)
    if 0 in (n, k):
        return Integer(0)
    n1 = n - 1

    # some special values
    if k == n1:
        return binomial(n, 2)
    if k == 2:
        return 2**n1 - 1

    # general recurrence
    return k*_stirling2(n1, k) + _stirling2(n1, k - 1)


def stirling(n, k, d=None, kind=2, signed=False):
    """Return Stirling number S(n, k) of the first or second (default) kind.

    The sum of all Stirling numbers of the second kind for k = 1
    through n is bell(n). The recurrence relationship for these numbers
    is::

    {0}       {n}   {0}      {n + 1}     {n}   {  n  }
    { } = 1;  { } = { } = 0; {     } = j*{ } + {     }
    {0}       {0}   {k}      {  k  }     {k}   {k - 1}

    where ``j`` is::
        ``n`` for Stirling numbers of the first kind
        ``-n`` for signed Stirling numbers of the first kind
        ``k`` for Stirling numbers of the second kind

    The first kind of Stirling number counts the number of permutations of
    ``n`` distinct items that have ``k`` cycles; the second kind counts the
    ways in which ``n`` distinct items can be partitioned into ``k`` parts.
    If ``d`` is given, the "reduced Stirling number of the second kind" is
    returned: ``S^{d}(n, k) = S(n - d + 1, k - d + 1)`` with ``n >= k >= d``.
    (This counts the ways to partition ``n`` consecutive integers into
    ``k`` groups with no pairwise difference less than ``d``. See example
    below.)

    To obtain the signed Stirling numbers of the first kind, use keyword
    ``signed=True``. Using this keyword automatically sets ``kind`` to 1.

    Examples
    ========

    >>> import itertools
    >>> from diofant.utilities.iterables import multiset_partitions

    First kind (unsigned by default):

    >>> [stirling(6, i, kind=1) for i in range(7)]
    [0, 120, 274, 225, 85, 15, 1]
    >>> perms = list(itertools.permutations(range(4)))
    >>> [sum(Permutation(p).cycles == i for p in perms) for i in range(5)]
    [0, 6, 11, 6, 1]
    >>> [stirling(4, i, kind=1) for i in range(5)]
    [0, 6, 11, 6, 1]

    First kind (signed):

    >>> [stirling(4, i, signed=True) for i in range(5)]
    [0, -6, 11, -6, 1]

    Second kind:

    >>> [stirling(10, i) for i in range(12)]
    [0, 1, 511, 9330, 34105, 42525, 22827, 5880, 750, 45, 1, 0]
    >>> sum(_) == bell(10)
    True
    >>> len(list(multiset_partitions(range(4), 2))) == stirling(4, 2)
    True

    Reduced second kind:

    >>> def delta(p):
    ...     if len(p) == 1:
    ...         return oo
    ...     return min(abs(i[0] - i[1]) for i in itertools.combinations(p, 2))
    >>> parts = multiset_partitions(range(5), 3)
    >>> d = 2
    >>> sum(1 for p in parts if all(delta(i) >= d for i in p))
    7
    >>> stirling(5, 3, 2)
    7

    References
    ==========

    * https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind
    * https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind

    See Also
    ========

    diofant.utilities.iterables.multiset_partitions

    """
    # TODO: make this a class like bell()

    n = as_int(n)
    k = as_int(k)
    if n < 0:
        raise ValueError('n must be nonnegative')
    if k > n:
        return Integer(0)
    if d:
        # assert k >= d
        # kind is ignored -- only kind=2 is supported
        return _stirling2(n - d + 1, k - d + 1)
    if signed:
        # kind is ignored -- only kind=1 is supported
        return (-1)**(n - k)*_stirling1(n, k)

    if kind == 1:
        return _stirling1(n, k)
    if kind == 2:
        return _stirling2(n, k)
    raise ValueError(f'kind must be 1 or 2, not {k}')