[docs]classfactorial(CombinatorialFunction):"""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 """deffdiff(self,argindex=1):fromsympyimportgamma,polygammaifargindex==1:returngamma(self.args[0]+1)*polygamma(0,self.args[0]+1)else:raiseArgumentIndexError(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=[]@classmethoddef_swing(cls,n):ifn<33:returncls._small_swing[n]else:N,primes=int(_sqrt(n)),[]forprimeinsieve.primerange(3,N+1):p,q=1,nwhileTrue:q//=primeifq>0:ifq&1==1:p*=primeelse:breakifp>1:primes.append(p)forprimeinsieve.primerange(N+1,n//3+1):if(n//prime)&1==1:primes.append(prime)L_product=R_product=1forprimeinsieve.primerange(n//2+1,n+1):L_product*=primeforprimeinprimes:R_product*=primereturnL_product*R_product@classmethoddef_recursive(cls,n):ifn<2:return1else:return(cls._recursive(n//2)**2)*cls._swing(n)@classmethoddefeval(cls,n):n=sympify(n)ifn.is_Number:ifnisS.Zero:returnS.OneelifnisS.Infinity:returnS.Infinityelifn.is_Integer:ifn.is_negative:returnS.ComplexInfinityelse:n=n.pifn<20:ifnotcls._small_factorials:result=1foriinrange(1,20):result*=icls._small_factorials.append(result)result=cls._small_factorials[n-1]# GMPY factorial is faster, use it when availableelifHAS_GMPY:fromsympy.core.compatibilityimportgmpyresult=gmpy.fac(n)else:bits=bin(n).count('1')result=cls._recursive(n)*2**(n-bits)returnInteger(result)def_eval_rewrite_as_gamma(self,n):fromsympyimportgammareturngamma(n+1)def_eval_rewrite_as_Product(self,n):fromsympyimportProductifn.is_nonnegativeandn.is_integer:i=Dummy('i',integer=True)returnProduct(i,(i,1,n))def_eval_is_integer(self):ifself.args[0].is_integerandself.args[0].is_nonnegative:returnTruedef_eval_is_positive(self):ifself.args[0].is_integerandself.args[0].is_nonnegative:returnTruedef_eval_is_even(self):x=self.args[0]ifx.is_integerandx.is_nonnegative:return(x-2).is_nonnegativedef_eval_is_composite(self):x=self.args[0]ifx.is_integerandx.is_nonnegative:return(x-3).is_nonnegativedef_eval_is_real(self):x=self.args[0]ifx.is_nonnegativeorx.is_noninteger:returnTruedef_eval_Mod(self,q):x=self.args[0]ifx.is_integerandx.is_nonnegativeandq.is_integer:aq=abs(q)d=x-aqifd.is_nonnegative:return0elifd==-1:''' Apply Wilson's theorem-if a natural number n > 1 is a prime number, (n-1)! = -1 mod n-and its inverse-if n > 4 is a composite number, (n-1)! = 0 mod n '''ifaq.is_prime:return-1%qelifaq.is_compositeand(aq-6).is_nonnegative:

[docs]classbinomial(CombinatorialFunction):"""Implementation of the binomial coefficient. It can be defined in two ways depending on its desired interpretation: C(n,k) = n!/(k!(n-k)!) or C(n, k) = 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/ """deffdiff(self,argindex=1):fromsympyimportpolygammaifargindex==1:# http://functions.wolfram.com/GammaBetaErf/Binomial/20/01/01/n,k=self.argsreturnbinomial(n,k)*(polygamma(0,n+1)- \
polygamma(0,n-k+1))elifargindex==2:# http://functions.wolfram.com/GammaBetaErf/Binomial/20/01/02/n,k=self.argsreturnbinomial(n,k)*(polygamma(0,n-k+1)- \
polygamma(0,k+1))else:raiseArgumentIndexError(self,argindex)@classmethoddef_eval(self,n,k):# n.is_Number and k.is_Integer and k != 1 and n != kfromsympy.functions.elementary.exponentialimportlogfromsympy.coreimportNifk.is_Integer:ifn.is_Integerandn>=0:n,k=int(n),int(k)ifk>n:returnS.Zeroelifk>n//2:k=n-kifHAS_GMPY:fromsympy.core.compatibilityimportgmpyreturnInteger(gmpy.bincoef(n,k))prime_count_estimate=N(n/log(n))# if the number of primes less than n is less than k, use prime sieve method# otherwise it is more memory efficient to compute factorials explicitlyifprime_count_estimate<k:M,result=int(_sqrt(n)),1forprimeinsieve.primerange(2,n+1):ifprime>n-k:result*=primeelifprime>n//2:continueelifprime>M:ifn%prime<k%prime:result*=primeelse:N,K=n,kexp=a=0whileN>0:a=int((N%prime)<(K%prime+a))N,K=N//prime,K//primeexp=a+expifexp>0:result*=prime**expelse:result=ff(n,k)/factorial(k)returnInteger(result)else:d=result=n-k+1foriinrange(2,k+1):d+=1result*=dresult/=ireturnresult@classmethoddefeval(cls,n,k):n,k=map(sympify,(n,k))d=n-kifd.is_zeroork.is_zero:returnS.Oneif(k-1).is_zero:returnnifk.is_integer:ifk.is_negative:returnS.Zeroifn.is_integerandn.is_nonnegativeandd.is_negative:returnS.Zeroifn.is_number:res=cls._eval(n,k)returnres.expand(basic=True)ifreselsereselifn.is_negativeandn.is_integerandnotk.is_integer:# a special case when binomial evaluates to complex infinityreturnS.ComplexInfinityelifk.is_number:fromsympyimportgammareturngamma(n+1)/(gamma(k+1)*gamma(n-k+1))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]ifn.is_Number:returnbinomial(*self.args)k=self.args[1]ifk.is_Addandnink.args:k=n-kifk.is_Integer:ifk==S.Zero:returnS.Oneelifk<0:returnS.Zeroelse:n=self.args[0]result=n-k+1foriinrange(2,k+1):result*=n-k+iresult/=ireturnresultelse:returnbinomial(*self.args)def_eval_rewrite_as_factorial(self,n,k):returnfactorial(n)/(factorial(k)*factorial(n-k))def_eval_rewrite_as_gamma(self,n,k):fromsympyimportgammareturngamma(n+1)/(gamma(k+1)*gamma(n-k+1))def_eval_rewrite_as_tractable(self,n,k):returnself._eval_rewrite_as_gamma(n,k).rewrite('tractable')def_eval_rewrite_as_FallingFactorial(self,n,k):ifk.is_integer:returnff(n,k)/factorial(k)def_eval_is_integer(self):n,k=self.argsifn.is_integerandk.is_integer:returnTrueelifk.is_integerisFalse:returnFalsedef_eval_is_nonnegative(self):ifself.args[0].is_integerandself.args[1].is_integer:ifself.args[0].is_nonnegativeandself.args[1].is_nonnegative: