Source code for sympy.functions.combinatorial.factorials

from __future__ import print_function, division

from sympy.core import S, sympify, Dummy, Mod
from sympy.core.function import Function, ArgumentIndexError
from sympy.core.logic import fuzzy_and
from sympy.core.numbers import Integer, pi
from sympy.core.relational import Eq

from sympy.ntheory import sieve

from math import sqrt as _sqrt

from sympy.core.compatibility import reduce, range, HAS_GMPY
from sympy.core.cache import cacheit

from sympy.polys.polytools import Poly

class CombinatorialFunction(Function):
    """Base class for combinatorial functions. """

    def _eval_simplify(self, ratio, measure, rational, inverse):
        from sympy.simplify.combsimp import combsimp
        # combinatorial function with non-integer arguments is
        # automatically passed to gammasimp
        expr = combsimp(self)
        if measure(expr) <= ratio*measure(self):
            return expr
        return self


###############################################################################
######################## FACTORIAL and MULTI-FACTORIAL ########################
###############################################################################


[docs]class factorial(CombinatorialFunction): r"""Implementation of factorial function over nonnegative integers. By convention (consistent with the gamma function and the binomial coefficients), factorial of a negative integer is complex infinity. The factorial is very important in combinatorics where it gives the number of ways in which `n` objects can be permuted. It also arises in calculus, probability, number theory, etc. There is strict relation of factorial with gamma function. In fact `n! = gamma(n+1)` for nonnegative integers. Rewrite of this kind is very useful in case of combinatorial simplification. Computation of the factorial is done using two algorithms. For small arguments a precomputed look up table is used. However for bigger input algorithm Prime-Swing is used. It is the fastest algorithm known and computes `n!` via prime factorization of special class of numbers, called here the 'Swing Numbers'. Examples ======== >>> from sympy import Symbol, factorial, S >>> n = Symbol('n', integer=True) >>> factorial(0) 1 >>> factorial(7) 5040 >>> factorial(-2) zoo >>> factorial(n) factorial(n) >>> factorial(2*n) factorial(2*n) >>> factorial(S(1)/2) factorial(1/2) See Also ======== factorial2, RisingFactorial, FallingFactorial """ def fdiff(self, argindex=1): from sympy import gamma, polygamma if argindex == 1: return gamma(self.args[0] + 1)*polygamma(0, self.args[0] + 1) else: raise ArgumentIndexError(self, argindex) _small_swing = [ 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003, 429, 6435, 6435, 109395, 12155, 230945, 46189, 969969, 88179, 2028117, 676039, 16900975, 1300075, 35102025, 5014575, 145422675, 9694845, 300540195, 300540195 ] _small_factorials = [] @classmethod def _swing(cls, n): if n < 33: return cls._small_swing[n] else: N, primes = int(_sqrt(n)), [] for prime in sieve.primerange(3, N + 1): p, q = 1, n while True: q //= prime if q > 0: if q & 1 == 1: p *= prime else: break if p > 1: primes.append(p) for prime in sieve.primerange(N + 1, n//3 + 1): if (n // prime) & 1 == 1: primes.append(prime) L_product = R_product = 1 for prime in sieve.primerange(n//2 + 1, n + 1): L_product *= prime for prime in primes: R_product *= prime return L_product*R_product @classmethod def _recursive(cls, n): if n < 2: return 1 else: return (cls._recursive(n//2)**2)*cls._swing(n) @classmethod def eval(cls, n): n = sympify(n) if n.is_Number: if n is S.Zero: return S.One elif n is S.Infinity: return S.Infinity elif n.is_Integer: if n.is_negative: return S.ComplexInfinity else: n = n.p if n < 20: if not cls._small_factorials: result = 1 for i in range(1, 20): result *= i cls._small_factorials.append(result) result = cls._small_factorials[n-1] # GMPY factorial is faster, use it when available elif HAS_GMPY: from sympy.core.compatibility import gmpy result = gmpy.fac(n) else: bits = bin(n).count('1') result = cls._recursive(n)*2**(n - bits) return Integer(result) def _facmod(self, n, q): res, N = 1, int(_sqrt(n)) # Exponent of prime p in n! is e_p(n) = [n/p] + [n/p**2] + ... # for p > sqrt(n), e_p(n) < sqrt(n), the primes with [n/p] = m, # occur consecutively and are grouped together in pw[m] for # simultaneous exponentiation at a later stage pw = [1]*N m = 2 # to initialize the if condition below for prime in sieve.primerange(2, n + 1): if m > 1: m, y = 0, n // prime while y: m += y y //= prime if m < N: pw[m] = pw[m]*prime % q else: res = res*pow(prime, m, q) % q for ex, bs in enumerate(pw): if ex == 0 or bs == 1: continue if bs == 0: return 0 res = res*pow(bs, ex, q) % q return res def _eval_Mod(self, q): n = self.args[0] if n.is_integer and n.is_nonnegative and q.is_integer: aq = abs(q) d = aq - n if d.is_nonpositive: return 0 else: isprime = aq.is_prime if d == 1: # Apply Wilson's theorem (if a natural number n > 1 # is a prime number, then (n-1)! = -1 mod n) and # its inverse (if n > 4 is a composite number, then # (n-1)! = 0 mod n) if isprime: return -1 % q elif isprime is False and (aq - 6).is_nonnegative: return 0 elif n.is_Integer and q.is_Integer: n, d, aq = map(int, (n, d, aq)) if isprime and (d - 1 < n): fc = self._facmod(d - 1, aq) fc = pow(fc, aq - 2, aq) if d%2: fc = -fc else: fc = self._facmod(n, aq) return Integer(fc % q) def _eval_rewrite_as_gamma(self, n, **kwargs): from sympy import gamma return gamma(n + 1) def _eval_rewrite_as_Product(self, n, **kwargs): from sympy import Product if n.is_nonnegative and n.is_integer: i = Dummy('i', integer=True) return Product(i, (i, 1, n)) def _eval_is_integer(self): if self.args[0].is_integer and self.args[0].is_nonnegative: return True def _eval_is_positive(self): if self.args[0].is_integer and self.args[0].is_nonnegative: return True def _eval_is_even(self): x = self.args[0] if x.is_integer and x.is_nonnegative: return (x - 2).is_nonnegative def _eval_is_composite(self): x = self.args[0] if x.is_integer and x.is_nonnegative: return (x - 3).is_nonnegative def _eval_is_real(self): x = self.args[0] if x.is_nonnegative or x.is_noninteger: return True
[docs]class MultiFactorial(CombinatorialFunction): pass
[docs]class subfactorial(CombinatorialFunction): r"""The subfactorial counts the derangements of n items and is defined for non-negative integers as: .. math:: !n = \begin{cases} 1 & n = 0 \\ 0 & n = 1 \\ (n-1)(!(n-1) + !(n-2)) & n > 1 \end{cases} It can also be written as ``int(round(n!/exp(1)))`` but the recursive definition with caching is implemented for this function. An interesting analytic expression is the following [2]_ .. math:: !x = \Gamma(x + 1, -1)/e which is valid for non-negative integers `x`. The above formula is not very useful incase of non-integers. :math:`\Gamma(x + 1, -1)` is single-valued only for integral arguments `x`, elsewhere on the positive real axis it has an infinite number of branches none of which are real. References ========== .. [1] https://en.wikipedia.org/wiki/Subfactorial .. [2] http://mathworld.wolfram.com/Subfactorial.html Examples ======== >>> from sympy import subfactorial >>> from sympy.abc import n >>> subfactorial(n + 1) subfactorial(n + 1) >>> subfactorial(5) 44 See Also ======== sympy.functions.combinatorial.factorials.factorial, sympy.utilities.iterables.generate_derangements, sympy.functions.special.gamma_functions.uppergamma """ @classmethod @cacheit def _eval(self, n): if not n: return S.One elif n == 1: return S.Zero return (n - 1)*(self._eval(n - 1) + self._eval(n - 2)) @classmethod def eval(cls, arg): if arg.is_Number: if arg.is_Integer and arg.is_nonnegative: return cls._eval(arg) elif arg is S.NaN: return S.NaN elif arg is S.Infinity: return S.Infinity def _eval_is_even(self): if self.args[0].is_odd and self.args[0].is_nonnegative: return True def _eval_is_integer(self): if self.args[0].is_integer and self.args[0].is_nonnegative: return True def _eval_rewrite_as_uppergamma(self, arg, **kwargs): from sympy import uppergamma return uppergamma(arg + 1, -1)/S.Exp1 def _eval_is_nonnegative(self): if self.args[0].is_integer and self.args[0].is_nonnegative: return True def _eval_is_odd(self): if self.args[0].is_even and self.args[0].is_nonnegative: return True
[docs]class factorial2(CombinatorialFunction): r"""The double factorial `n!!`, not to be confused with `(n!)!` The double factorial is defined for nonnegative integers and for odd negative integers as: .. math:: n!! = \begin{cases} 1 & n = 0 \\ n(n-2)(n-4) \cdots 1 & n\ \text{positive odd} \\ n(n-2)(n-4) \cdots 2 & n\ \text{positive even} \\ (n+2)!!/(n+2) & n\ \text{negative odd} \end{cases} References ========== .. [1] https://en.wikipedia.org/wiki/Double_factorial Examples ======== >>> from sympy import factorial2, var >>> var('n') n >>> factorial2(n + 1) factorial2(n + 1) >>> factorial2(5) 15 >>> factorial2(-1) 1 >>> factorial2(-5) 1/3 See Also ======== factorial, RisingFactorial, FallingFactorial """ @classmethod def eval(cls, arg): # TODO: extend this to complex numbers? if arg.is_Number: if not arg.is_Integer: raise ValueError("argument must be nonnegative integer " "or negative odd integer") # This implementation is faster than the recursive one # It also avoids "maximum recursion depth exceeded" runtime error if arg.is_nonnegative: if arg.is_even: k = arg / 2 return 2**k * factorial(k) return factorial(arg) / factorial2(arg - 1) if arg.is_odd: return arg*(S.NegativeOne)**((1 - arg)/2) / factorial2(-arg) raise ValueError("argument must be nonnegative integer " "or negative odd integer") def _eval_is_even(self): # Double factorial is even for every positive even input n = self.args[0] if n.is_integer: if n.is_odd: return False if n.is_even: if n.is_positive: return True if n.is_zero: return False def _eval_is_integer(self): # Double factorial is an integer for every nonnegative input, and for # -1 and -3 n = self.args[0] if n.is_integer: if (n + 1).is_nonnegative: return True if n.is_odd: return (n + 3).is_nonnegative def _eval_is_odd(self): # Double factorial is odd for every odd input not smaller than -3, and # for 0 n = self.args[0] if n.is_odd: return (n + 3).is_nonnegative if n.is_even: if n.is_positive: return False if n.is_zero: return True def _eval_is_positive(self): # Double factorial is positive for every nonnegative input, and for # every odd negative input which is of the form -1-4k for an # nonnegative integer k n = self.args[0] if n.is_integer: if (n + 1).is_nonnegative: return True if n.is_odd: return ((n + 1) / 2).is_even def _eval_rewrite_as_gamma(self, n, **kwargs): from sympy import gamma, Piecewise, sqrt return 2**(n/2)*gamma(n/2 + 1) * Piecewise((1, Eq(Mod(n, 2), 0)), (sqrt(2/pi), Eq(Mod(n, 2), 1)))
############################################################################### ######################## RISING and FALLING FACTORIALS ######################## ###############################################################################
[docs]class RisingFactorial(CombinatorialFunction): r""" Rising factorial (also called Pochhammer symbol) is a double valued function arising in concrete mathematics, hypergeometric functions and series expansions. It is defined by: .. math:: rf(x,k) = x \cdot (x+1) \cdots (x+k-1) where `x` can be arbitrary expression and `k` is an integer. For more information check "Concrete mathematics" by Graham, pp. 66 or visit http://mathworld.wolfram.com/RisingFactorial.html page. When `x` is a Poly instance of degree >= 1 with a single variable, `rf(x,k) = x(y) \cdot x(y+1) \cdots x(y+k-1)`, where `y` is the variable of `x`. This is as described in Peter Paule, "Greatest Factorial Factorization and Symbolic Summation", Journal of Symbolic Computation, vol. 20, pp. 235-268, 1995. Examples ======== >>> from sympy import rf, symbols, factorial, ff, binomial, Poly >>> from sympy.abc import x >>> n, k = symbols('n k', integer=True) >>> rf(x, 0) 1 >>> rf(1, 5) 120 >>> rf(x, 5) == x*(1 + x)*(2 + x)*(3 + x)*(4 + x) True >>> rf(Poly(x**3, x), 2) Poly(x**6 + 3*x**5 + 3*x**4 + x**3, x, domain='ZZ') Rewrite >>> rf(x, k).rewrite(ff) FallingFactorial(k + x - 1, k) >>> rf(x, k).rewrite(binomial) binomial(k + x - 1, k)*factorial(k) >>> rf(n, k).rewrite(factorial) factorial(k + n - 1)/factorial(n - 1) See Also ======== factorial, factorial2, FallingFactorial References ========== .. [1] https://en.wikipedia.org/wiki/Pochhammer_symbol """ @classmethod def eval(cls, x, k): x = sympify(x) k = sympify(k) if x is S.NaN or k is S.NaN: return S.NaN elif x is S.One: return factorial(k) elif k.is_Integer: if k is S.Zero: return S.One else: if k.is_positive: if x is S.Infinity: return S.Infinity elif x is S.NegativeInfinity: if k.is_odd: return S.NegativeInfinity else: return S.Infinity else: if isinstance(x, Poly): gens = x.gens if len(gens)!= 1: raise ValueError("rf only defined for " "polynomials on one generator") else: return reduce(lambda r, i: r*(x.shift(i).expand()), range(0, int(k)), 1) else: return reduce(lambda r, i: r*(x + i), range(0, int(k)), 1) else: if x is S.Infinity: return S.Infinity elif x is S.NegativeInfinity: return S.Infinity else: if isinstance(x, Poly): gens = x.gens if len(gens)!= 1: raise ValueError("rf only defined for " "polynomials on one generator") else: return 1/reduce(lambda r, i: r*(x.shift(-i).expand()), range(1, abs(int(k)) + 1), 1) else: return 1/reduce(lambda r, i: r*(x - i), range(1, abs(int(k)) + 1), 1) def _eval_rewrite_as_gamma(self, x, k, **kwargs): from sympy import gamma return gamma(x + k) / gamma(x) def _eval_rewrite_as_FallingFactorial(self, x, k, **kwargs): return FallingFactorial(x + k - 1, k) def _eval_rewrite_as_factorial(self, x, k, **kwargs): if x.is_integer and k.is_integer: return factorial(k + x - 1) / factorial(x - 1) def _eval_rewrite_as_binomial(self, x, k, **kwargs): if k.is_integer: return factorial(k) * binomial(x + k - 1, k) def _eval_is_integer(self): return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer, self.args[1].is_nonnegative)) def _sage_(self): import sage.all as sage return sage.rising_factorial(self.args[0]._sage_(), self.args[1]._sage_())
[docs]class FallingFactorial(CombinatorialFunction): r""" Falling factorial (related to rising factorial) is a double valued function arising in concrete mathematics, hypergeometric functions and series expansions. It is defined by .. math:: ff(x,k) = x \cdot (x-1) \cdots (x-k+1) where `x` can be arbitrary expression and `k` is an integer. For more information check "Concrete mathematics" by Graham, pp. 66 or visit http://mathworld.wolfram.com/FallingFactorial.html page. When `x` is a Poly instance of degree >= 1 with single variable, `ff(x,k) = x(y) \cdot x(y-1) \cdots x(y-k+1)`, where `y` is the variable of `x`. This is as described in Peter Paule, "Greatest Factorial Factorization and Symbolic Summation", Journal of Symbolic Computation, vol. 20, pp. 235-268, 1995. >>> from sympy import ff, factorial, rf, gamma, polygamma, binomial, symbols, Poly >>> from sympy.abc import x, k >>> n, m = symbols('n m', integer=True) >>> ff(x, 0) 1 >>> ff(5, 5) 120 >>> ff(x, 5) == x*(x-1)*(x-2)*(x-3)*(x-4) True >>> ff(Poly(x**2, x), 2) Poly(x**4 - 2*x**3 + x**2, x, domain='ZZ') >>> ff(n, n) factorial(n) Rewrite >>> ff(x, k).rewrite(gamma) (-1)**k*gamma(k - x)/gamma(-x) >>> ff(x, k).rewrite(rf) RisingFactorial(-k + x + 1, k) >>> ff(x, m).rewrite(binomial) binomial(x, m)*factorial(m) >>> ff(n, m).rewrite(factorial) factorial(n)/factorial(-m + n) See Also ======== factorial, factorial2, RisingFactorial References ========== .. [1] http://mathworld.wolfram.com/FallingFactorial.html """ @classmethod def eval(cls, x, k): x = sympify(x) k = sympify(k) if x is S.NaN or k is S.NaN: return S.NaN elif k.is_integer and x == k: return factorial(x) elif k.is_Integer: if k is S.Zero: return S.One else: if k.is_positive: if x is S.Infinity: return S.Infinity elif x is S.NegativeInfinity: if k.is_odd: return S.NegativeInfinity else: return S.Infinity else: if isinstance(x, Poly): gens = x.gens if len(gens)!= 1: raise ValueError("ff only defined for " "polynomials on one generator") else: return reduce(lambda r, i: r*(x.shift(-i).expand()), range(0, int(k)), 1) else: return reduce(lambda r, i: r*(x - i), range(0, int(k)), 1) else: if x is S.Infinity: return S.Infinity elif x is S.NegativeInfinity: return S.Infinity else: if isinstance(x, Poly): gens = x.gens if len(gens)!= 1: raise ValueError("rf only defined for " "polynomials on one generator") else: return 1/reduce(lambda r, i: r*(x.shift(i).expand()), range(1, abs(int(k)) + 1), 1) else: return 1/reduce(lambda r, i: r*(x + i), range(1, abs(int(k)) + 1), 1) def _eval_rewrite_as_gamma(self, x, k, **kwargs): from sympy import gamma return (-1)**k*gamma(k - x) / gamma(-x) def _eval_rewrite_as_RisingFactorial(self, x, k, **kwargs): return rf(x - k + 1, k) def _eval_rewrite_as_binomial(self, x, k, **kwargs): if k.is_integer: return factorial(k) * binomial(x, k) def _eval_rewrite_as_factorial(self, x, k, **kwargs): if x.is_integer and k.is_integer: return factorial(x) / factorial(x - k) def _eval_is_integer(self): return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer, self.args[1].is_nonnegative)) def _sage_(self): import sage.all as sage return sage.falling_factorial(self.args[0]._sage_(), self.args[1]._sage_())
rf = RisingFactorial ff = FallingFactorial ############################################################################### ########################### BINOMIAL COEFFICIENTS ############################# ###############################################################################
[docs]class binomial(CombinatorialFunction): r"""Implementation of the binomial coefficient. It can be defined in two ways depending on its desired interpretation: .. math:: \binom{n}{k} = \frac{n!}{k!(n-k)!}\ \text{or}\ \binom{n}{k} = \frac{ff(n, k)}{k!} First, in a strict combinatorial sense it defines the number of ways we can choose `k` elements from a set of `n` elements. In this case both arguments are nonnegative integers and binomial is computed using an efficient algorithm based on prime factorization. The other definition is generalization for arbitrary `n`, however `k` must also be nonnegative. This case is very useful when evaluating summations. For the sake of convenience for negative integer `k` this function will return zero no matter what valued is the other argument. To expand the binomial when `n` is a symbol, use either ``expand_func()`` or ``expand(func=True)``. The former will keep the polynomial in factored form while the latter will expand the polynomial itself. See examples for details. Examples ======== >>> from sympy import Symbol, Rational, binomial, expand_func >>> n = Symbol('n', integer=True, positive=True) >>> binomial(15, 8) 6435 >>> binomial(n, -1) 0 Rows of Pascal's triangle can be generated with the binomial function: >>> for N in range(8): ... print([binomial(N, i) for i in range(N + 1)]) ... [1] [1, 1] [1, 2, 1] [1, 3, 3, 1] [1, 4, 6, 4, 1] [1, 5, 10, 10, 5, 1] [1, 6, 15, 20, 15, 6, 1] [1, 7, 21, 35, 35, 21, 7, 1] As can a given diagonal, e.g. the 4th diagonal: >>> N = -4 >>> [binomial(N, i) for i in range(1 - N)] [1, -4, 10, -20, 35] >>> binomial(Rational(5, 4), 3) -5/128 >>> binomial(Rational(-5, 4), 3) -195/128 >>> binomial(n, 3) binomial(n, 3) >>> binomial(n, 3).expand(func=True) n**3/6 - n**2/2 + n/3 >>> expand_func(binomial(n, 3)) n*(n - 2)*(n - 1)/6 References ========== .. [1] https://www.johndcook.com/blog/binomial_coefficients/ """ def fdiff(self, argindex=1): from sympy import polygamma if argindex == 1: # http://functions.wolfram.com/GammaBetaErf/Binomial/20/01/01/ n, k = self.args return binomial(n, k)*(polygamma(0, n + 1) - \ polygamma(0, n - k + 1)) elif argindex == 2: # http://functions.wolfram.com/GammaBetaErf/Binomial/20/01/02/ n, k = self.args return binomial(n, k)*(polygamma(0, n - k + 1) - \ polygamma(0, k + 1)) else: raise ArgumentIndexError(self, argindex) @classmethod def _eval(self, n, k): # n.is_Number and k.is_Integer and k != 1 and n != k from sympy.functions.elementary.exponential import log from sympy.core import N if k.is_Integer: if n.is_Integer and n >= 0: n, k = int(n), int(k) if k > n: return S.Zero elif k > n // 2: k = n - k if HAS_GMPY: from sympy.core.compatibility import gmpy return Integer(gmpy.bincoef(n, k)) d, result = n - k, 1 for i in range(1, k + 1): d += 1 result = result * d // i return Integer(result) else: d, result = n - k, 1 for i in range(1, k + 1): d += 1 result *= d result /= i return result @classmethod def eval(cls, n, k): n, k = map(sympify, (n, k)) d = n - k n_nonneg, n_isint = n.is_nonnegative, n.is_integer if k.is_zero or ((n_nonneg or n_isint is False) and d.is_zero): return S.One if (k - 1).is_zero or ((n_nonneg or n_isint is False) and (d - 1).is_zero): return n if k.is_integer: if k.is_negative or (n_nonneg and n_isint and d.is_negative): return S.Zero elif n.is_number: res = cls._eval(n, k) return res.expand(basic=True) if res else res elif n_nonneg is False and n_isint: # a special case when binomial evaluates to complex infinity return S.ComplexInfinity elif k.is_number: from sympy import gamma return gamma(n + 1)/(gamma(k + 1)*gamma(n - k + 1)) def _eval_Mod(self, q): n, k = self.args if any(x.is_integer is False for x in (n, k, q)): raise ValueError("Integers expected for binomial Mod") if all(x.is_Integer for x in (n, k, q)): n, k = map(int, (n, k)) aq, res = abs(q), 1 # handle negative integers k or n if k < 0: return 0 if n < 0: n = -n + k - 1 res = -1 if k%2 else 1 # non negative integers k and n if k > n: return 0 isprime = aq.is_prime aq = int(aq) if isprime: if aq < n: # use Lucas Theorem N, K = n, k while N or K: res = res*binomial(N % aq, K % aq) % aq N, K = N // aq, K // aq else: # use Factorial Modulo d = n - k if k > d: k, d = d, k kf = 1 for i in range(2, k + 1): kf = kf*i % aq df = kf for i in range(k + 1, d + 1): df = df*i % aq res *= df for i in range(d + 1, n + 1): res = res*i % aq res *= pow(kf*df % aq, aq - 2, aq) res %= aq else: # Binomial Factorization is performed by calculating the # exponents of primes <= n in `n! /(k! (n - k)!)`, # for non-negative integers n and k. As the exponent of # prime in n! is e_p(n) = [n/p] + [n/p**2] + ... # the exponent of prime in binomial(n, k) would be # e_p(n) - e_p(k) - e_p(n - k) M = int(_sqrt(n)) for prime in sieve.primerange(2, n + 1): if prime > n - k: res = res*prime % aq elif prime > n // 2: continue elif prime > M: if n % prime < k % prime: res = res*prime % aq else: N, K = n, k exp = a = 0 while N > 0: a = int((N % prime) < (K % prime + a)) N, K = N // prime, K // prime exp += a if exp > 0: res *= pow(prime, exp, aq) res %= aq return Integer(res % q) def _eval_expand_func(self, **hints): """ Function to expand binomial(n, k) when m is positive integer Also, n is self.args[0] and k is self.args[1] while using binomial(n, k) """ n = self.args[0] if n.is_Number: return binomial(*self.args) k = self.args[1] if k.is_Add and n in k.args: k = n - k if k.is_Integer: if k == S.Zero: return S.One elif k < 0: return S.Zero else: n, result = self.args[0], 1 for i in range(1, k + 1): result *= n - k + i result /= i return result else: return binomial(*self.args) def _eval_rewrite_as_factorial(self, n, k, **kwargs): return factorial(n)/(factorial(k)*factorial(n - k)) def _eval_rewrite_as_gamma(self, n, k, **kwargs): from sympy import gamma return gamma(n + 1)/(gamma(k + 1)*gamma(n - k + 1)) def _eval_rewrite_as_tractable(self, n, k, **kwargs): return self._eval_rewrite_as_gamma(n, k).rewrite('tractable') def _eval_rewrite_as_FallingFactorial(self, n, k, **kwargs): if k.is_integer: return ff(n, k) / factorial(k) def _eval_is_integer(self): n, k = self.args if n.is_integer and k.is_integer: return True elif k.is_integer is False: return False def _eval_is_nonnegative(self): n, k = self.args if n.is_integer and k.is_integer: if n.is_nonnegative or k.is_negative or k.is_even: return True elif k.is_even is False: return False