Navigation

Source code for sympy.core.numbers

fromcoreimportCfromsympifyimportconverter,sympify,_sympify,SympifyErrorfrombasicimportBasicfromsingletonimportS,SingletonfromexprimportExpr,AtomicExprfromdecoratorsimport_sympifyit,deprecatedfromcacheimportcacheit,clear_cacheimportsympy.mpmathasmpmathimportsympy.mpmath.libmpasmlibfromsympy.mpmath.libmpimportmpf_pow,mpf_pi,mpf_e,phi_fixedfromsympy.mpmath.ctx_mpimportmpnumericimportdecimalrnd=mlib.round_nearest# TODO: we should use the warnings module_errdict={"divide":False}

_gcdcache={}# TODO caching with decorator, but not to degrade performance

[docs]defigcd(a,b):"""Computes positive, integer greatest common divisor of two numbers. The algorithm is based on the well known Euclid's algorithm. To improve speed, igcd() has its own caching mechanism implemented. """try:return_gcdcache[(a,b)]exceptKeyError:ifaandb:ifb<0:b=-bwhileb:a,b=b,a%belse:a=abs(aorb)_gcdcache[(a,b)]=areturna

[docs]defilcm(a,b):"""Computes integer least common multiple of two numbers. """ifa==0andb==0:return0else:returna*b//igcd(a,b)

[docs]classNumber(AtomicExpr):""" Represents any kind of number in sympy. Floating point numbers are represented by the Float class. Integer numbers (of any size), together with rational numbers (again, there is no limit on their size) are represented by the Rational class. If you want to represent, for example, ``1+sqrt(2)``, then you need to do:: Rational(1) + sqrt(Rational(2)) """is_commutative=Trueis_comparable=Trueis_bounded=Trueis_finite=Trueis_number=True__slots__=[]# Used to make max(x._prec, y._prec) return x._prec when only x is a float_prec=-1is_Number=Truedef__new__(cls,*obj):iflen(obj)==1:obj=obj[0]ifisinstance(obj,(int,long)):returnInteger(obj)ifisinstance(obj,tuple)andlen(obj)==2:returnRational(*obj)ifisinstance(obj,(float,mpmath.mpf,decimal.Decimal)):returnFloat(obj)ifisinstance(obj,str):val=sympify(obj)ifisinstance(val,Number):returnvalelse:raiseValueError('String "%s" does not denote a Number'%obj)ifisinstance(obj,Number):returnobjraiseTypeError("expected str|int|long|float|Decimal|Number object but got %r"%(obj))def_as_mpf_val(self,prec):"""Evaluation of mpf tuple accurate to at least prec bits."""raiseNotImplementedError('%s needs ._as_mpf_val() method'% \
(self.__class__.__name__))def_eval_evalf(self,prec):returnFloat._new(self._as_mpf_val(prec),prec)def_as_mpf_op(self,prec):prec=max(prec,self._prec)returnself._as_mpf_val(prec),precdef__float__(self):returnmlib.to_float(self._as_mpf_val(53))def_eval_conjugate(self):returnselfdef_eval_order(self,*symbols):# Order(5, x, y) -> Order(1,x,y)returnC.Order(S.One,*symbols)@classmethoddefclass_key(cls):return1,0,'Number'defsort_key(self,order=None):returnself.class_key(),(0,()),(),selfdef__eq__(self,other):raiseNotImplementedError('%s needs .__eq__() method'%(self.__class__.__name__))def__ne__(self,other):raiseNotImplementedError('%s needs .__ne__() method'%(self.__class__.__name__))def__lt__(self,other):raiseNotImplementedError('%s needs .__lt__() method'%(self.__class__.__name__))def__le__(self,other):raiseNotImplementedError('%s needs .__le__() method'%(self.__class__.__name__))def__gt__(self,other):return_sympify(other).__lt__(self)def__ge__(self,other):return_sympify(other).__le__(self)def__hash__(self):returnsuper(Number,self).__hash__()@propertydefis_number(self):returnTruedefas_coeff_mul(self,*deps):# a -> c * tifself.is_Rational:returnself,tuple()elifself.is_negative:returnS.NegativeOne,(-self,)returnS.One,(self,)defas_coeff_add(self,*deps):# a -> c + tifself.is_Rational:returnself,tuple()returnS.Zero,(self,)

[docs]deflimit_denominator(self,max_denominator=1000000):"""Closest Rational to self with denominator at most max_denominator. >>> from sympy import Rational >>> Rational('3.141592653589793').limit_denominator(10) 22/7 >>> Rational('3.141592653589793').limit_denominator(100) 311/99 """# Algorithm notes: For any real number x, define a *best upper# approximation* to x to be a rational number p/q such that:## (1) p/q >= x, and# (2) if p/q > r/s >= x then s > q, for any rational r/s.## Define *best lower approximation* similarly. Then it can be# proved that a rational number is a best upper or lower# approximation to x if, and only if, it is a convergent or# semiconvergent of the (unique shortest) continued fraction# associated to x.## To find a best rational approximation with denominator <= M,# we find the best upper and lower approximations with# denominator <= M and take whichever of these is closer to x.# In the event of a tie, the bound with smaller denominator is# chosen. If both denominators are equal (which can happen# only when max_denominator == 1 and self is midway between# two integers) the lower bound---i.e., the floor of self, is# taken.ifmax_denominator<1:raiseValueError("max_denominator should be at least 1")ifself.q<=max_denominator:returnselfp0,q0,p1,q1=0,1,1,0n,d=self.p,self.qwhileTrue:a=n//dq2=q0+a*q1ifq2>max_denominator:breakp0,q0,p1,q1=p1,q1,p0+a*p1,q2n,d=d,n-a*dk=(max_denominator-q0)//q1bound1=Rational(p0+k*p1,q0+k*q1)bound2=Rational(p1,q1)ifabs(bound2-self)<=abs(bound1-self):returnbound2else:returnbound1

[docs]deffactors(self,limit=None,use_trial=True,use_rho=False,use_pm1=False,verbose=False):"""A wrapper to factorint which return factors of self that are smaller than limit (or cheap to compute). Special methods of factoring are disabled by default so that only trial division is used. """fromsympy.ntheoryimportfactorintf=factorint(self.p,limit=limit,use_trial=use_trial,use_rho=use_rho,use_pm1=use_pm1,verbose=verbose).copy()forp,einfactorint(self.q,limit=limit,use_trial=use_trial,use_rho=use_rho,use_pm1=use_pm1,verbose=verbose).items():try:f[p]+=-eexceptKeyError:f[p]=-eiflen(f)>1and1inf:delf[1]returnf

[docs]classInteger(Rational):q=1is_integer=Trueis_Integer=True__slots__=['p']def_as_mpf_val(self,prec):returnmlib.from_int(self.p)def_mpmath_(self,prec,rnd):returnmpmath.make_mpf(self._as_mpf_val(prec))# TODO caching with decorator, but not to degrade performance@int_tracedef__new__(cls,i):ival=int(i)try:return_intcache[ival]exceptKeyError:# We only work with well-behaved integer types. This converts, for# example, numpy.int32 instances.ifival==0:obj=S.Zeroelifival==1:obj=S.Oneelifival==-1:obj=S.NegativeOneelse:obj=Expr.__new__(cls)obj.p=ival_intcache[ival]=objreturnobjdef__getnewargs__(self):return(self.p,)# Arithmetic operations are here for efficiencydef__int__(self):returnself.pdef__neg__(self):returnInteger(-self.p)def__abs__(self):ifself.p>=0:returnselfelse:returnInteger(-self.p)def__divmod__(self,other):returndivmod(self.p,other.p)# TODO make it decorator + bytecodehacks?def__add__(a,b):ifisinstance(b,(int,long)):returnInteger(a.p+b)elifisinstance(b,Integer):returnInteger(a.p+b.p)returnRational.__add__(a,b)# a,b -not- b,adef__radd__(a,b):ifisinstance(b,(int,long)):returnInteger(b+a.p)elifisinstance(b,Integer):returnInteger(b.p+a.p)returnRational.__add__(a,b)def__sub__(a,b):ifisinstance(b,(int,long)):returnInteger(a.p-b)elifisinstance(b,Integer):returnInteger(a.p-b.p)returnRational.__sub__(a,b)def__rsub__(a,b):ifisinstance(b,(int,long)):returnInteger(b-a.p)elifisinstance(b,Integer):returnInteger(b.p-a.p)returnRational.__rsub__(a,b)def__mul__(a,b):ifisinstance(b,(int,long)):returnInteger(a.p*b)elifisinstance(b,Integer):returnInteger(a.p*b.p)returnRational.__mul__(a,b)def__rmul__(a,b):ifisinstance(b,(int,long)):returnInteger(b*a.p)elifisinstance(b,Integer):returnInteger(b.p*a.p)returnRational.__mul__(a,b)def__mod__(a,b):ifisinstance(b,(int,long)):returnInteger(a.p%b)elifisinstance(b,Integer):returnInteger(a.p%b.p)returnRational.__mod__(a,b)def__rmod__(a,b):ifisinstance(b,(int,long)):returnInteger(b%a.p)elifisinstance(b,Integer):returnInteger(b.p%a.p)returnRational.__rmod__(a,b)def__eq__(a,b):ifisinstance(b,(int,long)):return(a.p==b)elifisinstance(b,Integer):return(a.p==b.p)returnRational.__eq__(a,b)def__ne__(a,b):ifisinstance(b,(int,long)):return(a.p!=b)elifisinstance(b,Integer):return(a.p!=b.p)returnRational.__ne__(a,b)def__gt__(a,b):ifisinstance(b,(int,long)):return(a.p>b)elifisinstance(b,Integer):return(a.p>b.p)returnRational.__gt__(a,b)def__lt__(a,b):ifisinstance(b,(int,long)):return(a.p<b)elifisinstance(b,Integer):return(a.p<b.p)returnRational.__lt__(a,b)def__ge__(a,b):ifisinstance(b,(int,long)):return(a.p>=b)elifisinstance(b,Integer):return(a.p>=b.p)returnRational.__ge__(a,b)def__le__(a,b):ifisinstance(b,(int,long)):return(a.p<=b)elifisinstance(b,Integer):return(a.p<=b.p)returnRational.__le__(a,b)def__hash__(self):returnsuper(Integer,self).__hash__()def__index__(self):returnself.p########################################def_eval_is_odd(self):returnbool(self.p%2)def_eval_power(b,e):""" Tries to do some simplifications on b ** e, where b is an instance of Integer Returns None if no further simplifications can be done When exponent is a fraction (so we have for example a square root), we try to find a simpler representation by factoring the argument up to factors of 2**15, e.g. - 4**Rational(1,2) becomes 2 - (-4)**Rational(1,2) becomes 2*I - (2**(3+7)*3**(6+7))**Rational(1,7) becomes 6*18**(3/7) Further simplification would require a special call to factorint on the argument which is not done here for sake of speed. """fromsympyimportperfect_powerifeisS.NaN:returnS.NaNifbisS.One:returnS.OneifbisS.NegativeOne:returnifeisS.Infinity:ifb>S.One:returnS.InfinityifbisS.NegativeOne:returnS.NaN# cases for 0 and 1 are done in their respective classesreturnS.Infinity+S.ImaginaryUnit*S.Infinityifnotisinstance(e,Number):# simplify when exp is even# (-2) ** k --> 2 ** kc,t=b.as_coeff_mul()ife.is_evenandisinstance(c,Number)andc<0:return(-c*Mul(*t))**eifnotisinstance(e,Rational):returnifeisS.Halfandb<0:# we extract I for this special case since everyone is doing soreturnS.ImaginaryUnit*Pow(-b,e)ife<0:# invert base and change sign on exponentne=-eifb<0:ife.q!=1:return-(S.NegativeOne)**((e.p%e.q)/S(e.q))*Rational(1,-b)**neelse:return(S.NegativeOne)**ne*Rational(1,-b)**neelse:returnRational(1,b)**ne# see if base is a perfect root, sqrt(4) --> 2b_pos=int(abs(b))x,xexact=integer_nthroot(b_pos,e.q)ifxexact:# if it's a perfect root we've finishedresult=Integer(x**abs(e.p))ifb<0:result*=(-1)**ereturnresult# The following is an algorithm where we collect perfect roots# from the factors of base.# if it's not an nth root, it still might be a perfect powerp=perfect_power(b_pos)ifpisnotFalse:dict={p[0]:p[1]}else:dict=Integer(b_pos).factors(limit=2**15)# now process the dict of factorsifb.is_negative:dict[-1]=1out_int=1# integer partout_rad=1# extracted radicalssqr_int=1sqr_gcd=0sqr_dict={}forprime,exponentindict.items():exponent*=e.p# remove multiples of e.q, e.g. (2**12)**(1/10) -> 2*(2**2)**(1/10)div_e,div_m=divmod(exponent,e.q)ifdiv_e>0:out_int*=prime**div_eifdiv_m>0:# see if the reduced exponent shares a gcd with e.q# (2**2)**(1/10) -> 2**(1/5)g=igcd(div_m,e.q)ifg!=1:out_rad*=Pow(prime,Rational(div_m//g,e.q//g))else:sqr_dict[prime]=div_m# identify gcd of remaining powersforp,exinsqr_dict.iteritems():ifsqr_gcd==0:sqr_gcd=exelse:sqr_gcd=igcd(sqr_gcd,ex)ifsqr_gcd==1:breakfork,vinsqr_dict.iteritems():sqr_int*=k**(v//sqr_gcd)ifsqr_int==bandout_int==1andout_rad==1:result=Noneelse:result=out_int*out_rad*Pow(sqr_int,Rational(sqr_gcd,e.q))returnresultdef_eval_is_prime(self):ifself.p<0:returnFalsedefas_numer_denom(self):returnself,S.Onedef__floordiv__(self,other):returnInteger(self.p//Integer(other).p)def__rfloordiv__(self,other):returnInteger(Integer(other).p//self.p)