[docs]@timethis('combsimp')defcombsimp(expr):r""" Simplify combinatorial expressions. This function takes as input an expression containing factorials, binomials, Pochhammer symbol and other "combinatorial" functions, and tries to minimize the number of those functions and reduce the size of their arguments. The algorithm works by rewriting all combinatorial functions as expressions involving rising factorials (Pochhammer symbols) and applies recurrence relations and other transformations applicable to rising factorials, to reduce their arguments, possibly letting the resulting rising factorial to cancel. Rising factorials with the second argument being an integer are expanded into polynomial forms and finally all other rising factorial are rewritten in terms of more familiar functions. If the initial expression consisted of gamma functions alone, the result is expressed in terms of gamma functions. If the initial expression consists of gamma function with some other combinatorial, the result is expressed in terms of gamma functions. If the result is expressed using gamma functions, the following three additional steps are performed: 1. Reduce the number of gammas by applying the reflection theorem gamma(x)*gamma(1-x) == pi/sin(pi*x). 2. Reduce the number of gammas by applying the multiplication theorem gamma(x)*gamma(x+1/n)*...*gamma(x+(n-1)/n) == C*gamma(n*x). 3. Reduce the number of prefactors by absorbing them into gammas, where possible. All transformation rules can be found (or was derived from) here: 1. http://functions.wolfram.com/GammaBetaErf/Pochhammer/17/01/02/ 2. http://functions.wolfram.com/GammaBetaErf/Pochhammer/27/01/0005/ Examples ======== >>> from sympy.simplify import combsimp >>> from sympy import factorial, binomial >>> from sympy.abc import n, k >>> combsimp(factorial(n)/factorial(n - 3)) n*(n - 2)*(n - 1) >>> combsimp(binomial(n+1, k+1)/binomial(n, k)) (n + 1)/(k + 1) """# as a rule of thumb, if the expression contained gammas initially, it# probably makes sense to retain themas_gamma=expr.has(gamma)as_factorial=expr.has(factorial)as_binomial=expr.has(binomial)expr=expr.replace(binomial,lambdan,k:_rf((n-k+1).expand(),k.expand())/_rf(1,k.expand()))expr=expr.replace(factorial,lambdan:_rf(1,n.expand()))expr=expr.rewrite(gamma)expr=expr.replace(gamma,lambdan:_rf(1,(n-1).expand()))ifas_gamma:expr=expr.replace(_rf,lambdaa,b:gamma(a+b)/gamma(a))else:expr=expr.replace(_rf,lambdaa,b:binomial(a+b-1,b)*gamma(b+1))defrule(n,k):coeff,rewrite=S.One,Falsecn,_n=n.as_coeff_Add()if_nandcn.is_Integerandcn:coeff*=_rf(_n+1,cn)/_rf(_n-k+1,cn)rewrite=Truen=_n# this sort of binomial has already been removed by# rising factorials but is left here in case the order# of rule application is changedifk.is_Add:ck,_k=k.as_coeff_Add()if_kandck.is_Integerandck:coeff*=_rf(n-ck-_k+1,ck)/_rf(_k+1,ck)rewrite=Truek=_kifrewrite:returncoeff*binomial(n,k)expr=expr.replace(binomial,rule)defrule_gamma(expr,level=0):""" Simplify products of gamma functions further. """ifexpr.is_Atom:returnexprdefgamma_rat(x):# helper to simplify ratios of gammaswas=x.count(gamma)xx=x.replace(gamma,lambdan:_rf(1,(n-1).expand()).replace(_rf,lambdaa,b:gamma(a+b)/gamma(a)))ifxx.count(gamma)<was:x=xxreturnxdefgamma_factor(x):# return True if there is a gamma factor in shallow argsifx.funcisgamma:returnTrueifx.is_Addorx.is_Mul:returnany(gamma_factor(xi)forxiinx.args)ifx.is_Powand(x.exp.is_integerorx.base.is_positive):returngamma_factor(x.base)returnFalse# recursion stepiflevel==0:expr=expr.func(*[rule_gamma(x,level+1)forxinexpr.args])level+=1ifnotexpr.is_Mul:returnexpr# non-commutative stepiflevel==1:args,nc=expr.args_cnc()ifnotargs:returnexprifnc:returnrule_gamma(Mul._from_args(args),level+1)*Mul._from_args(nc)level+=1# pure gamma handling, not factor absorbtioniflevel==2:sifted=sift(expr.args,gamma_factor)gamma_ind=Mul(*sifted.pop(False,[]))d=Mul(*sifted.pop(True,[]))assertnotsiftednd,dd=d.as_numer_denom()foripassinrange(2):args=list(ordered(Mul.make_args(nd)))fori,niinenumerate(args):ifni.is_Add:ni,dd=Add(*[rule_gamma(gamma_rat(a/dd),level+1)forainni.args]).as_numer_denom()args[i]=niifnotdd.has(gamma):breaknd=Mul(*args)ifipass==0andnotgamma_factor(nd):breaknd,dd=dd,nd# now process in reversed orderexpr=gamma_ind*nd/ddifnot(expr.is_Muland(gamma_factor(dd)orgamma_factor(nd))):returnexprlevel+=1# iteration until constantiflevel==3:whileTrue:was=exprexpr=rule_gamma(expr,4)ifexpr==was:returnexprnumer_gammas=[]denom_gammas=[]numer_others=[]denom_others=[]defexplicate(p):ifpisS.One:returnNone,[]b,e=p.as_base_exp()ife.is_Integer:ifb.funcisgamma:returnTrue,[b.args[0]]*eelse:returnFalse,[b]*eelse:returnFalse,[p]newargs=list(ordered(expr.args))whilenewargs:n,d=newargs.pop().as_numer_denom()isg,l=explicate(n)ifisg:numer_gammas.extend(l)elifisgisFalse:numer_others.extend(l)isg,l=explicate(d)ifisg:denom_gammas.extend(l)elifisgisFalse:denom_others.extend(l)# =========== level 2 work: pure gamma manipulation =========# Try to reduce the number of gamma factors by applying the# reflection formula gamma(x)*gamma(1-x) = pi/sin(pi*x)forgammas,numer,denomin[(numer_gammas,numer_others,denom_others),(denom_gammas,denom_others,numer_others)]:new=[]whilegammas:g1=gammas.pop()ifg1.is_integer:new.append(g1)continuefori,g2inenumerate(gammas):n=g1+g2-1ifnotn.is_Integer:continuenumer.append(S.Pi)denom.append(sin(S.Pi*g1))gammas.pop(i)ifn>0:forkinrange(n):numer.append(1-g1+k)elifn<0:forkinrange(-n):denom.append(-g1-k)breakelse:new.append(g1)# /!\ updating IN PLACEgammas[:]=new# Try to reduce the number of gammas by using the duplication# theorem to cancel an upper and lower: gamma(2*s)/gamma(s) =# 2**(2*s + 1)/(4*sqrt(pi))*gamma(s + 1/2). Although this could# be done with higher argument ratios like gamma(3*x)/gamma(x),# this would not reduce the number of gammas as in this case.forng,dg,no,doin[(numer_gammas,denom_gammas,numer_others,denom_others),(denom_gammas,numer_gammas,denom_others,numer_others)]:whileTrue:forxinng:foryindg:n=x-2*yifn.is_Integer:breakelse:continuebreakelse:breakng.remove(x)dg.remove(y)ifn>0:forkinrange(n):no.append(2*y+k)elifn<0:forkinrange(-n):do.append(2*y-1-k)ng.append(y+S(1)/2)no.append(2**(2*y-1))do.append(sqrt(S.Pi))# Try to reduce the number of gamma factors by applying the# multiplication theorem (used when n gammas with args differing# by 1/n mod 1 are encountered).## run of 2 with args differing by 1/2## >>> combsimp(gamma(x)*gamma(x+S.Half))# 2*sqrt(2)*2**(-2*x - 1/2)*sqrt(pi)*gamma(2*x)## run of 3 args differing by 1/3 (mod 1)## >>> combsimp(gamma(x)*gamma(x+S(1)/3)*gamma(x+S(2)/3))# 6*3**(-3*x - 1/2)*pi*gamma(3*x)# >>> combsimp(gamma(x)*gamma(x+S(1)/3)*gamma(x+S(5)/3))# 2*3**(-3*x - 1/2)*pi*(3*x + 2)*gamma(3*x)#def_run(coeffs):# find runs in coeffs such that the difference in terms (mod 1)# of t1, t2, ..., tn is 1/nu=list(uniq(coeffs))foriinrange(len(u)):dj=([((u[j]-u[i])%1,j)forjinrange(i+1,len(u))])forone,jindj:ifone.p==1andone.q!=1:n=one.qgot=[i]get=list(range(1,n))ford,jindj:m=n*difm.is_Integerandminget:get.remove(m)got.append(j)ifnotget:breakelse:continuefori,jinenumerate(got):c=u[j]coeffs.remove(c)got[i]=creturnone.q,got[0],got[1:]def_mult_thm(gammas,numer,denom):# pull off and analyze the leading coefficient from each gamma arg# looking for runs in those Rationals# expr -> coeff + resid -> rats[resid] = coeffrats={}forgingammas:c,resid=g.as_coeff_Add()rats.setdefault(resid,[]).append(c)# look for runs in Rationals for each residkeys=sorted(rats,key=default_sort_key)forresidinkeys:coeffs=list(sorted(rats[resid]))new=[]whileTrue:run=_run(coeffs)ifrunisNone:break# process the sequence that was found:# 1) convert all the gamma functions to have the right# argument (could be off by an integer)# 2) append the factors corresponding to the theorem# 3) append the new gamma functionn,ui,other=run# (1)foruinother:con=resid+u-1forkinrange(int(u-ui)):numer.append(con-k)con=n*(resid+ui)# for (2) and (3)# (2)numer.append((2*S.Pi)**(S(n-1)/2)*n**(S(1)/2-con))# (3)new.append(con)# restore resid to coeffsrats[resid]=[resid+cforcincoeffs]+new# rebuild the gamma argumentsg=[]forresidinkeys:g+=rats[resid]# /!\ updating IN PLACEgammas[:]=gforl,numer,denomin[(numer_gammas,numer_others,denom_others),(denom_gammas,denom_others,numer_others)]:_mult_thm(l,numer,denom)# =========== level >= 2 work: factor absorbtion =========iflevel>=2:# Try to absorb factors into the gammas: x*gamma(x) -> gamma(x + 1)# and gamma(x)/(x - 1) -> gamma(x - 1)# This code (in particular repeated calls to find_fuzzy) can be very# slow.deffind_fuzzy(l,x):ifnotl:returnS1,T1=compute_ST(x)foryinl:S2,T2=inv[y]ifT1!=T2or(notS1.intersection(S2)and(S1!=set()orS2!=set())):continue# XXX we want some simplification (e.g. cancel or# simplify) but no matter what it's slow.a=len(cancel(x/y).free_symbols)b=len(x.free_symbols)c=len(y.free_symbols)# TODO is there a better heuristic?ifa==0and(b>0orc>0):returny# We thus try to avoid expensive calls by building the following# "invariants": For every factor or gamma function argument# - the set of free symbols S# - the set of functional components T# We will only try to absorb if T1==T2 and (S1 intersect S2 != emptyset# or S1 == S2 == emptyset)inv={}defcompute_ST(expr):ifexprininv:returninv[expr]return(expr.free_symbols,expr.atoms(Function).union(set(e.expforeinexpr.atoms(Pow))))defupdate_ST(expr):inv[expr]=compute_ST(expr)forexprinnumer_gammas+denom_gammas+numer_others+denom_others:update_ST(expr)forgammas,numer,denomin[(numer_gammas,numer_others,denom_others),(denom_gammas,denom_others,numer_others)]:new=[]whilegammas:g=gammas.pop()cont=Truewhilecont:cont=Falsey=find_fuzzy(numer,g)ifyisnotNone:numer.remove(y)ify!=g:numer.append(y/g)update_ST(y/g)g+=1cont=Truey=find_fuzzy(denom,g-1)ifyisnotNone:denom.remove(y)ify!=g-1:numer.append((g-1)/y)update_ST((g-1)/y)g-=1cont=Truenew.append(g)# /!\ updating IN PLACEgammas[:]=new# =========== rebuild expr ==================================returnMul(*[gamma(g)forginnumer_gammas]) \
/Mul(*[gamma(g)forgindenom_gammas]) \
*Mul(*numer_others)/Mul(*denom_others)# (for some reason we cannot use Basic.replace in this case)was=factor(expr)expr=rule_gamma(was)ifexpr!=was:expr=factor(expr)ifnotas_gamma:ifas_factorial:expr=expr.rewrite(factorial)elifas_binomial:expr=expr.rewrite(binomial)returnexpr