Navigation

Source code for sympy.series.gruntz

"""Limits======Implemented according to the PhD thesishttp://www.cybertester.com/data/gruntz.pdf, which contains very thoroughdescriptions of the algorithm including many examples. We summarize herethe gist of it.All functions are sorted according to how rapidly varying they are atinfinity using the following rules. Any two functions f and g can becompared using the properties of L:L=lim log|f(x)| / log|g(x)| (for x -> oo)We define >, < ~ according to:: 1. f > g .... L=+-oo we say that: - f is greater than any power of g - f is more rapidly varying than g - f goes to infinity/zero faster than g 2. f < g .... L=0 we say that: - f is lower than any power of g 3. f ~ g .... L!=0, +-oo we say that: - both f and g are bounded from above and below by suitable integral powers of the otherExamples========:: 2 < x < exp(x) < exp(x**2) < exp(exp(x)) 2 ~ 3 ~ -5 x ~ x**2 ~ x**3 ~ 1/x ~ x**m ~ -x exp(x) ~ exp(-x) ~ exp(2x) ~ exp(x)**2 ~ exp(x+exp(-x)) f ~ 1/fSo we can divide all the functions into comparability classes (x and x^2belong to one class, exp(x) and exp(-x) belong to some other class). Inprinciple, we could compare any two functions, but in our algorithm, wedon't compare anything below the class 2~3~-5 (for example log(x) isbelow this), so we set 2~3~-5 as the lowest comparability class.Given the function f, we find the list of most rapidly varying (mrv set)subexpressions of it. This list belongs to the same comparability class.Let's say it is {exp(x), exp(2x)}. Using the rule f ~ 1/f we find anelement "w" (either from the list or a new one) from the samecomparability class which goes to zero at infinity. In our example weset w=exp(-x) (but we could also set w=exp(-2x) or w=exp(-3x) ...). Werewrite the mrv set using w, in our case {1/w, 1/w^2}, and substitute itinto f. Then we expand f into a series in w:: f = c0*w^e0 + c1*w^e1 + ... + O(w^en), where e0<e1<...<en, c0!=0but for x->oo, lim f = lim c0*w^e0, because all the other terms go to zero,because w goes to zero faster than the ci and ei. So:: for e0>0, lim f = 0 for e0<0, lim f = +-oo (the sign depends on the sign of c0) for e0=0, lim f = lim c0We need to recursively compute limits at several places of the algorithm, butas is shown in the PhD thesis, it always finishes.Important functions from the implementation:compare(a, b, x) compares "a" and "b" by computing the limit L.mrv(e, x) returns list of most rapidly varying (mrv) subexpressions of "e"rewrite(e, Omega, x, wsym) rewrites "e" in terms of wleadterm(f, x) returns the lowest power term in the series of fmrv_leadterm(e, x) returns the lead term (c0, e0) for elimitinf(e, x) computes lim e (for x->oo)limit(e, z, z0) computes any limit by converting it to the case x->ooAll the functions are really simple and straightforward exceptrewrite(), which is the most difficult/complex part of the algorithm.When the algorithm fails, the bugs are usually in the series expansion(i.e. in SymPy) or in rewrite.This code is almost exact rewrite of the Maple code inside the Gruntzthesis.Debugging---------Because the gruntz algorithm is highly recursive, it's difficult tofigure out what went wrong inside a debugger. Instead, turn on nicedebug prints by defining the environment variable SYMPY_DEBUG. Forexample:[user@localhost]: SYMPY_DEBUG=True ./bin/isympyIn [1]: limit(sin(x)/x, x, 0)limitinf(_x*sin(1/_x), _x) = 1+-mrv_leadterm(_x*sin(1/_x), _x) = (1, 0)| +-mrv(_x*sin(1/_x), _x) = set([_x])| | +-mrv(_x, _x) = set([_x])| | +-mrv(sin(1/_x), _x) = set([_x])| | +-mrv(1/_x, _x) = set([_x])| | +-mrv(_x, _x) = set([_x])| +-mrv_leadterm(exp(_x)*sin(exp(-_x)), _x, set([exp(_x)])) = (1, 0)| +-rewrite(exp(_x)*sin(exp(-_x)), set([exp(_x)]), _x, _w) = (1/_w*sin(_w), -_x)| +-sign(_x, _x) = 1| +-mrv_leadterm(1, _x) = (1, 0)+-sign(0, _x) = 0+-limitinf(1, _x) = 1And check manually which line is wrong. Then go to the source code anddebug this function to figure out the exact problem."""from__future__importprint_function,divisionfromsympy.coreimportBasic,S,oo,Symbol,I,Dummy,Wild,Mulfromsympy.functionsimportlog,expfromsympy.series.orderimportOrderfromsympy.simplifyimportpowsimpfromsympyimportcacheitfromsympy.core.compatibilityimportreducefromsympy.utilities.timeutilsimporttimethistimeit=timethis('gruntz')fromsympy.utilities.miscimportdebug_decoratorasdebug

[docs]defcompare(a,b,x):"""Returns "<" if a<b, "=" for a == b, ">" for a>b"""# log(exp(...)) must always be simplified here for terminationla,lb=log(a),log(b)ifisinstance(a,Basic)anda.funcisexp:la=a.args[0]ifisinstance(b,Basic)andb.funcisexp:lb=b.args[0]c=limitinf(la/lb,x)ifc==0:return"<"elifc.is_infinite:return">"else:return"="

[docs]classSubsSet(dict):""" Stores (expr, dummy) pairs, and how to rewrite expr-s. The gruntz algorithm needs to rewrite certain expressions in term of a new variable w. We cannot use subs, because it is just too smart for us. For example:: > Omega=[exp(exp(_p - exp(-_p))/(1 - 1/_p)), exp(exp(_p))] > O2=[exp(-exp(_p) + exp(-exp(-_p))*exp(_p)/(1 - 1/_p))/_w, 1/_w] > e = exp(exp(_p - exp(-_p))/(1 - 1/_p)) - exp(exp(_p)) > e.subs(Omega[0],O2[0]).subs(Omega[1],O2[1]) -1/w + exp(exp(p)*exp(-exp(-p))/(1 - 1/p)) is really not what we want! So we do it the hard way and keep track of all the things we potentially want to substitute by dummy variables. Consider the expression:: exp(x - exp(-x)) + exp(x) + x. The mrv set is {exp(x), exp(-x), exp(x - exp(-x))}. We introduce corresponding dummy variables d1, d2, d3 and rewrite:: d3 + d1 + x. This class first of all keeps track of the mapping expr->variable, i.e. will at this stage be a dictionary:: {exp(x): d1, exp(-x): d2, exp(x - exp(-x)): d3}. [It turns out to be more convenient this way round.] But sometimes expressions in the mrv set have other expressions from the mrv set as subexpressions, and we need to keep track of that as well. In this case, d3 is really exp(x - d2), so rewrites at this stage is:: {d3: exp(x-d2)}. The function rewrite uses all this information to correctly rewrite our expression in terms of w. In this case w can be choosen to be exp(-x), i.e. d2. The correct rewriting then is:: exp(-w)/w + 1/w + x. """def__init__(self):self.rewrites={}def__repr__(self):returnsuper(SubsSet,self).__repr__()+', '+self.rewrites.__repr__()def__getitem__(self,key):ifnotkeyinself:self[key]=Dummy()returndict.__getitem__(self,key)defdo_subs(self,e):forexpr,varinself.items():e=e.subs(var,expr)returne

[docs]defmeets(self,s2):"""Tell whether or not self and s2 have non-empty intersection"""returnset(self.keys()).intersection(list(s2.keys()))!=set()

[docs]defunion(self,s2,exps=None):"""Compute the union of self and s2, adjusting exps"""res=self.copy()tr={}forexpr,varins2.items():ifexprinself:ifexps:exps=exps.subs(var,res[expr])tr[var]=res[expr]else:res[expr]=varforvar,rewrins2.rewrites.items():res.rewrites[var]=rewr.subs(tr)returnres,exps

[docs]defmrv(e,x):"""Returns a SubsSet of most rapidly varying (mrv) subexpressions of 'e', and e rewritten in terms of these"""e=powsimp(e,deep=True,combine='exp')ifnotisinstance(e,Basic):raiseTypeError("e should be an instance of Basic")ifnote.has(x):returnSubsSet(),eelife==x:s=SubsSet()returns,s[x]elife.is_Mulore.is_Add:i,d=e.as_independent(x)# throw away x-independent termsifd.func!=e.func:s,expr=mrv(d,x)returns,e.func(i,expr)a,b=d.as_two_terms()s1,e1=mrv(a,x)s2,e2=mrv(b,x)returnmrv_max1(s1,s2,e.func(i,e1,e2),x)elife.is_Pow:b,e=e.as_base_exp()ife.has(x):returnmrv(exp(e*log(b)),x)else:s,expr=mrv(b,x)returns,expr**eelife.funcislog:s,expr=mrv(e.args[0],x)returns,log(expr)elife.funcisexp:# We know from the theory of this algorithm that exp(log(...)) may always# be simplified here, and doing so is vital for termination.ife.args[0].funcislog:returnmrv(e.args[0].args[0],x)# if a product has an infinite factor the result will be# infinite if there is no zero, otherwise NaN; here, we# consider the result infinite if any factor is infiniteli=limitinf(e.args[0],x)ifany(_.is_infinitefor_inMul.make_args(li)):s1=SubsSet()e1=s1[e]s2,e2=mrv(e.args[0],x)su=s1.union(s2)[0]su.rewrites[e1]=exp(e2)returnmrv_max3(s1,e1,s2,exp(e2),su,e1,x)else:s,expr=mrv(e.args[0],x)returns,exp(expr)elife.is_Function:l=[mrv(a,x)foraine.args]l2=[sfor(s,_)inlifs!=SubsSet()]iflen(l2)!=1:# e.g. something like BesselJ(x, x)raiseNotImplementedError("MRV set computation for functions in"" several variables not implemented.")s,ss=l2[0],SubsSet()args=[ss.do_subs(x[1])forxinl]returns,e.func(*args)elife.is_Derivative:raiseNotImplementedError("MRV set computation for derviatives"" not implemented yet.")returnmrv(e.args[0],x)raiseNotImplementedError("Don't know how to calculate the mrv of '%s'"%e)

[docs]defmrv_max3(f,expsf,g,expsg,union,expsboth,x):"""Computes the maximum of two sets of expressions f and g, which are in the same comparability class, i.e. max() compares (two elements of) f and g and returns either (f, expsf) [if f is larger], (g, expsg) [if g is larger] or (union, expsboth) [if f, g are of the same class]. """ifnotisinstance(f,SubsSet):raiseTypeError("f should be an instance of SubsSet")ifnotisinstance(g,SubsSet):raiseTypeError("g should be an instance of SubsSet")iff==SubsSet():returng,expsgelifg==SubsSet():returnf,expsfeliff.meets(g):returnunion,expsbothc=compare(list(f.keys())[0],list(g.keys())[0],x)ifc==">":returnf,expsfelifc=="<":returng,expsgelse:ifc!="=":raiseValueError("c should be =")returnunion,expsboth

[docs]defmrv_max1(f,g,exps,x):"""Computes the maximum of two sets of expressions f and g, which are in the same comparability class, i.e. mrv_max1() compares (two elements of) f and g and returns the set, which is in the higher comparability class of the union of both, if they have the same order of variation. Also returns exps, with the appropriate substitutions made. """u,b=f.union(g,exps)returnmrv_max3(f,g.do_subs(exps),g,f.do_subs(exps),u,b,x)

@debug@cacheit@timeit

[docs]defsign(e,x):""" Returns a sign of an expression e(x) for x->oo. :: e > 0 for x sufficiently large ... 1 e == 0 for x sufficiently large ... 0 e < 0 for x sufficiently large ... -1 The result of this function is currently undefined if e changes sign arbitarily often for arbitrarily large x (e.g. sin(x)). Note that this returns zero only if e is *constantly* zero for x sufficiently large. [If e is constant, of course, this is just the same thing as the sign of e.] """fromsympyimportsignas_signifnotisinstance(e,Basic):raiseTypeError("e should be an instance of Basic")ife.is_positive:return1elife.is_negative:return-1elife.is_zero:return0elifnote.has(x):return_sign(e)elife==x:return1elife.is_Mul:a,b=e.as_two_terms()sa=sign(a,x)ifnotsa:return0returnsa*sign(b,x)elife.funcisexp:return1elife.is_Pow:s=sign(e.base,x)ifs==1:return1ife.exp.is_Integer:returns**e.expelife.funcislog:returnsign(e.args[0]-1,x)# if all else fails, do it the hard wayc0,e0=mrv_leadterm(e,x)returnsign(c0,x)

@debug@timeit@cacheit

[docs]deflimitinf(e,x):"""Limit e(x) for x-> oo"""#rewrite e in terms of tractable functions onlye=e.rewrite('tractable',deep=True)ifnote.has(x):returne# e is a constantife.has(Order):e=e.expand().removeO()ifnotx.is_positive:# We make sure that x.is_positive is True so we# get all the correct mathematical behavior from the expression.# We need a fresh variable.p=Dummy('p',positive=True,finite=True)e=e.subs(x,p)x=pc0,e0=mrv_leadterm(e,x)sig=sign(e0,x)ifsig==1:returnS.Zero# e0>0: lim f = 0elifsig==-1:# e0<0: lim f = +-oo (the sign depends on the sign of c0)ifc0.match(I*Wild("a",exclude=[I])):returnc0*oos=sign(c0,x)#the leading term shouldn't be 0:ifs==0:raiseValueError("Leading term should not be 0")returns*ooelifsig==0:returnlimitinf(c0,x)# e0=0: lim f = lim c0

[docs]defcalculate_series(e,x,logx=None):""" Calculates at least one term of the series of "e" in "x". This is a place that fails most often, so it is in its own function. """fromsympy.polysimportcancelfortine.lseries(x,logx=logx):t=cancel(t)ift.simplify():breakreturnt

@debug@timeit@cacheit

[docs]defmrv_leadterm(e,x):"""Returns (c0, e0) for e."""Omega=SubsSet()ifnote.has(x):return(e,S.Zero)ifOmega==SubsSet():Omega,exps=mrv(e,x)ifnotOmega:# e really does not depend on x after simplificationseries=calculate_series(e,x)c0,e0=series.leadterm(x)ife0!=0:raiseValueError("e0 should be 0")returnc0,e0ifxinOmega:#move the whole omega up (exponentiate each term):Omega_up=moveup2(Omega,x)e_up=moveup([e],x)[0]exps_up=moveup([exps],x)[0]# NOTE: there is no need to move this down!e=e_upOmega=Omega_upexps=exps_up## The positive dummy, w, is used here so log(w*2) etc. will expand;# a unique dummy is needed in this algorithm## For limits of complex functions, the algorithm would have to be# improved, or just find limits of Re and Im components separately.#w=Dummy("w",real=True,positive=True,finite=True)f,logw=rewrite(exps,Omega,x,w)series=calculate_series(f,w,logx=logw)returnseries.leadterm(w)

[docs]defbuild_expression_tree(Omega,rewrites):r""" Helper function for rewrite. We need to sort Omega (mrv set) so that we replace an expression before we replace any expression in terms of which it has to be rewritten:: e1 ---> e2 ---> e3 \ -> e4 Here we can do e1, e2, e3, e4 or e1, e2, e4, e3. To do this we assemble the nodes into a tree, and sort them by height. This function builds the tree, rewrites then sorts the nodes. """classNode:defht(self):returnreduce(lambdax,y:x+y,[x.ht()forxinself.before],1)nodes={}forexpr,vinOmega:n=Node()n.before=[]n.var=vn.expr=exprnodes[v]=nfor_,vinOmega:ifvinrewrites:n=nodes[v]r=rewrites[v]for_,v2inOmega:ifr.has(v2):n.before.append(nodes[v2])returnnodes

@debug@timeit

[docs]defrewrite(e,Omega,x,wsym):"""e(x) ... the function Omega ... the mrv set wsym ... the symbol which is going to be used for w Returns the rewritten e in terms of w and log(w). See test_rewrite1() for examples and correct results. """fromsympyimportilcmifnotisinstance(Omega,SubsSet):raiseTypeError("Omega should be an instance of SubsSet")iflen(Omega)==0:raiseValueError("Length can not be 0")#all items in Omega must be exponentialsfortinOmega.keys():ifnott.funcisexp:raiseValueError("Value should be exp")rewrites=Omega.rewritesOmega=list(Omega.items())nodes=build_expression_tree(Omega,rewrites)Omega.sort(key=lambdax:nodes[x[1]].ht(),reverse=True)# make sure we know the sign of each exp() term; after the loop,# g is going to be the "w" - the simplest one in the mrv setforg,_inOmega:sig=sign(g.args[0],x)ifsig!=1andsig!=-1:raiseNotImplementedError('Result depends on the sign of %s'%sig)ifsig==1:wsym=1/wsym# if g goes to oo, substitute 1/w#O2 is a list, which results by rewriting each item in Omega using "w"O2=[]denominators=[]forf,varinOmega:c=limitinf(f.args[0]/g.args[0],x)ifc.is_Rational:denominators.append(c.q)arg=f.args[0]ifvarinrewrites:ifnotrewrites[var].funcisexp:raiseValueError("Value should be exp")arg=rewrites[var].args[0]O2.append((var,exp((arg-c*g.args[0]).expand())*wsym**c))#Remember that Omega contains subexpressions of "e". So now we find#them in "e" and substitute them for our rewriting, stored in O2# the following powsimp is necessary to automatically combine exponentials,# so that the .subs() below succeeds:# TODO this should not be necessaryf=powsimp(e,deep=True,combine='exp')fora,binO2:f=f.subs(a,b)for_,varinOmega:assertnotf.has(var)#finally compute the logarithm of w (logw).logw=g.args[0]ifsig==1:logw=-logw# log(w)->log(1/w)=-log(w)# Some parts of sympy have difficulty computing series expansions with# non-integral exponents. The following heuristic improves the situation:exponent=reduce(ilcm,denominators,1)f=f.subs(wsym,wsym**exponent)logw/=exponentreturnf,logw

[docs]defgruntz(e,z,z0,dir="+"):""" Compute the limit of e(z) at the point z0 using the Gruntz algorithm. z0 can be any expression, including oo and -oo. For dir="+" (default) it calculates the limit from the right (z->z0+) and for dir="-" the limit from the left (z->z0-). For infinite z0 (oo or -oo), the dir argument doesn't matter. This algorithm is fully described in the module docstring in the gruntz.py file. It relies heavily on the series expansion. Most frequently, gruntz() is only used if the faster limit() function (which uses heuristics) fails. """ifnotisinstance(z,Symbol):raiseNotImplementedError("Second argument must be a Symbol")#convert all limits to the limit z->oo; sign of z is handled in limitinfr=Noneifz0==oo:r=limitinf(e,z)elifz0==-oo:r=limitinf(e.subs(z,-z),z)else:ifstr(dir)=="-":e0=e.subs(z,z0-1/z)elifstr(dir)=="+":e0=e.subs(z,z0+1/z)else:raiseNotImplementedError("dir must be '+' or '-'")r=limitinf(e0,z)# This is a bit of a heuristic for nice results... we always rewrite# tractable functions in terms of familiar intractable ones.# It might be nicer to rewrite the exactly to what they were initially,# but that would take some work to implement.returnr.rewrite('intractable',deep=True)