Source code for sympy.core.numbers

from core import C
from sympify import converter, sympify, _sympify, SympifyError
from basic import Basic
from singleton import S, Singleton
from expr import Expr, AtomicExpr
from decorators import _sympifyit, deprecated
from cache import cacheit, clear_cache
import sympy.mpmath as mpmath
import sympy.mpmath.libmp as mlib
from sympy.mpmath.libmp import mpf_pow, mpf_pi, mpf_e, phi_fixed
from sympy.mpmath.ctx_mp import mpnumeric

import decimal


rnd = mlib.round_nearest


# TODO: we should use the warnings module
_errdict = {"divide": False}
[docs]def seterr(divide=False): """ Should sympy raise an exception on 0/0 or return a nan? divide == True .... raise an exception divide == False ... return nan """ if _errdict["divide"] != divide: clear_cache() _errdict["divide"] = divide # (a,b) -> gcd(a,b)
_gcdcache = {} # TODO caching with decorator, but not to degrade performance
[docs]def igcd(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)] except KeyError: if a and b: if b < 0: b = -b while b: a, b = b, a % b else: a = abs(a or b) _gcdcache[(a,b)] = a return a
[docs]def ilcm(a, b): """Computes integer least common multiple of two numbers. """ if a == 0 and b == 0: return 0 else: return a * b // igcd(a, b)
def igcdex(a, b): """Returns x, y, g such that g = x*a + y*b = gcd(a, b). >>> from sympy.core.numbers import igcdex >>> igcdex(2, 3) (-1, 1, 1) >>> igcdex(10, 12) (-1, 1, 2) >>> x, y, g = igcdex(100, 2004) >>> x, y, g (-20, 1, 4) >>> x*100 + y*2004 4 """ if (not a) and (not b): return (0, 1, 0) if not a: return (0, b//abs(b), abs(b)) if not b: return (a//abs(a), 0, abs(a)) if a < 0: a, x_sign = -a, -1 else: x_sign = 1 if b < 0: b, y_sign = -b, -1 else: y_sign = 1 x, y, r, s = 1, 0, 0, 1 while b: (c, q) = (a % b, a // b) (a, b, r, s, x, y) = (b, c, x-q*r, y-q*s, r, s) return (x*x_sign, y*y_sign, a)
[docs]class Number(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 = True is_comparable = True is_bounded = True is_finite = True is_number = True __slots__ = [] # Used to make max(x._prec, y._prec) return x._prec when only x is a float _prec = -1 is_Number = True def __new__(cls, *obj): if len(obj)==1: obj=obj[0] if isinstance(obj, (int, long)): return Integer(obj) if isinstance(obj, tuple) and len(obj) == 2: return Rational(*obj) if isinstance(obj, (float, mpmath.mpf, decimal.Decimal)): return Float(obj) if isinstance(obj, str): val = sympify(obj) if isinstance(val, Number): return val else: raise ValueError('String "%s" does not denote a Number'%obj) if isinstance(obj, Number): return obj raise TypeError("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.""" raise NotImplementedError('%s needs ._as_mpf_val() method' % \ (self.__class__.__name__)) def _eval_evalf(self, prec): return Float._new(self._as_mpf_val(prec), prec) def _as_mpf_op(self, prec): prec = max(prec, self._prec) return self._as_mpf_val(prec), prec def __float__(self): return mlib.to_float(self._as_mpf_val(53)) def _eval_conjugate(self): return self def _eval_order(self, *symbols): # Order(5, x, y) -> Order(1,x,y) return C.Order(S.One, *symbols) @classmethod def class_key(cls): return 1, 0, 'Number' def sort_key(self, order=None): return self.class_key(), (0, ()), (), self def __eq__(self, other): raise NotImplementedError('%s needs .__eq__() method' % (self.__class__.__name__)) def __ne__(self, other): raise NotImplementedError('%s needs .__ne__() method' % (self.__class__.__name__)) def __lt__(self, other): raise NotImplementedError('%s needs .__lt__() method' % (self.__class__.__name__)) def __le__(self, other): raise NotImplementedError('%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): return super(Number, self).__hash__() @property def is_number(self): return True def as_coeff_mul(self, *deps): # a -> c * t if self.is_Rational: return self, tuple() elif self.is_negative: return S.NegativeOne, (-self,) return S.One, (self,) def as_coeff_add(self, *deps): # a -> c + t if self.is_Rational: return self, tuple() return S.Zero, (self,)
[docs] def gcd(self, other): """Compute greatest common divisor of input arguments. """ _ = _sympify(other) return S.One
[docs] def lcm(self, other): """Compute least common multiple of input arguments. """ other = _sympify(other) return self*other
[docs] def cofactors(self, other): """Compute GCD and cofactors of input arguments. """ other = _sympify(other) return S.One, self, other
[docs] def as_coeff_Mul(self): """Efficiently extract the coefficient of a product. """ return self, S.One
[docs]class Float(Number): """ Represents a floating point number. It is capable of representing arbitrary-precision floating-point numbers **Usage** :: Float(3.5) 3.5 # (the 3.5 was converted from a python float) Float("3.0000000000000005") >>> from sympy import Float >>> Float((1,3,0,2)) # mpmath tuple: (-1)**1 * 3 * 2**0; 3 has 2 bits -3.00000000000000 **Notes** - Float(x) with x being a Python int/long will return Integer(x) """ is_real = True is_irrational = False is_integer = False __slots__ = ['_mpf_', '_prec'] # mpz can't be pickled def __getnewargs__(self): return (mlib.to_pickable(self._mpf_),) def __getstate__(self): d = Expr.__getstate__(self).copy() del d["_mpf_"] return mlib.to_pickable(self._mpf_), d def __setstate__(self, state): _mpf_, d = state _mpf_ = mlib.from_pickable(_mpf_) self._mpf_ = _mpf_ Expr.__setstate__(self, d) is_Float = True def floor(self): return C.Integer(int(mlib.to_int(mlib.mpf_floor(self._mpf_, self._prec)))) def ceiling(self): return C.Integer(int(mlib.to_int(mlib.mpf_ceil(self._mpf_, self._prec)))) @property def num(self): return mpmath.mpf(self._mpf_) def _as_mpf_val(self, prec): return self._mpf_ def _as_mpf_op(self, prec): return self._mpf_, max(prec, self._prec) def __new__(cls, num, prec=15): prec = mlib.libmpf.dps_to_prec(prec) if isinstance(num, (int, long)): return Integer(num) if isinstance(num, (str, decimal.Decimal)): _mpf_ = mlib.from_str(str(num), prec, rnd) elif isinstance(num, tuple) and len(num) == 4: if type(num[1]) is str: # it's a hexadecimal (coming from a pickled object) # assume that it is in standard form num = list(num) num[1] = long(num[1], 16) _mpf_ = tuple(num) else: _mpf_ = mpmath.mpf( S.NegativeOne ** num[0] * num[1] * 2 ** num[2])._mpf_ else: _mpf_ = mpmath.mpf(num)._mpf_ if not num: return C.Zero() obj = Expr.__new__(cls) obj._mpf_ = _mpf_ obj._prec = prec return obj @classmethod def _new(cls, _mpf_, _prec): if _mpf_ == mlib.fzero: return S.Zero obj = Expr.__new__(cls) obj._mpf_ = _mpf_ obj._prec = _prec return obj def _hashable_content(self): return (self._mpf_, self._prec) def _eval_is_positive(self): return self.num > 0 def _eval_is_negative(self): return self.num < 0 def __neg__(self): return Float._new(mlib.mpf_neg(self._mpf_), self._prec) @_sympifyit('other', NotImplemented) def __mul__(self, other): if isinstance(other, Number): rhs, prec = other._as_mpf_op(self._prec) return Float._new(mlib.mpf_mul(self._mpf_, rhs, prec, rnd), prec) return Number.__mul__(self, other) @_sympifyit('other', NotImplemented) def __mod__(self, other): if isinstance(other, Number): rhs, prec = other._as_mpf_op(self._prec) return Float._new(mlib.mpf_mod(self._mpf_, rhs, prec, rnd), prec) return Number.__mod__(self, other) @_sympifyit('other', NotImplemented) def __rmod__(self, other): if isinstance(other, Number): rhs, prec = other._as_mpf_op(self._prec) return Float._new(mlib.mpf_mod(rhs, self._mpf_, prec, rnd), prec) return Number.__rmod__(self, other) @_sympifyit('other', NotImplemented) def __add__(self, other): if (other is S.NaN) or (self is NaN): return S.NaN if isinstance(other, Number): rhs, prec = other._as_mpf_op(self._prec) return Float._new(mlib.mpf_add(self._mpf_, rhs, prec, rnd), prec) return Number.__add__(self, other) def _eval_power(self, e): """ e is symbolic object but not equal to 0, 1 (-p) ** r -> exp(r * log(-p)) -> exp(r * (log(p) + I*Pi)) -> -> p ** r * (sin(Pi*r) + cos(Pi*r) * I) """ if isinstance(e, Number): if isinstance(e, Integer): prec = self._prec return Float._new(mlib.mpf_pow_int(self._mpf_, e.p, prec, rnd), prec) e, prec = e._as_mpf_op(self._prec) b = self._mpf_ try: y = mpf_pow(b, e, prec, rnd) return Float._new(y, prec) except mlib.ComplexResult: re, im = mlib.mpc_pow((b, mlib.fzero), (e, mlib.fzero), prec, rnd) return Float._new(re, prec) + Float._new(im, prec) * S.ImaginaryUnit def __abs__(self): return Float._new(mlib.mpf_abs(self._mpf_), self._prec) def __int__(self): return int(mlib.to_int(self._mpf_)) def __eq__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy != other --> not == if isinstance(other, NumberSymbol): if other.is_irrational: return False return other.__eq__(self) if isinstance(other, FunctionClass): #cos as opposed to cos(x) return False if isinstance(other, Number): return bool(mlib.mpf_eq(self._mpf_, other._as_mpf_val(self._prec))) return False # Float != non-Number def __ne__(self, other): try: other = _sympify(other) except SympifyError: return True # sympy != other if isinstance(other, NumberSymbol): if other.is_irrational: return True return other.__ne__(self) if isinstance(other, FunctionClass): #cos as opposed to cos(x) return True if isinstance(other, Number): return bool(not mlib.mpf_eq(self._mpf_, other._as_mpf_val(self._prec))) return True # Float != non-Number def __lt__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy > other if isinstance(other, NumberSymbol): return other.__ge__(self) if other.is_comparable: other = other.evalf() if isinstance(other, Number): return bool(mlib.mpf_lt(self._mpf_, other._as_mpf_val(self._prec))) return Expr.__lt__(self, other) def __le__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy > other --> ! <= if isinstance(other, NumberSymbol): return other.__gt__(self) if other.is_comparable: other = other.evalf() if isinstance(other, Number): return bool(mlib.mpf_le(self._mpf_, other._as_mpf_val(self._prec))) return Expr.__le__(self, other) def __hash__(self): return super(Float, self).__hash__() def epsilon_eq(self, other, epsilon="10e-16"): return abs(self - other) < Float(epsilon) def _sage_(self): import sage.all as sage return sage.RealNumber(str(self)) # Add sympify converters
converter[float] = converter[decimal.Decimal] = Float # this is here to work nicely in Sage RealNumber = Float @deprecated
[docs]def Real(*args, **kwargs): # pragma: no cover """Deprecated alias for the Float constructor.""" return Float(*args, **kwargs)
[docs]class Rational(Number): """Represents integers and rational numbers (p/q) of any size. **Examples** >>> from sympy import Rational >>> from sympy.abc import x, y >>> Rational(3) 3 >>> Rational(1,2) 1/2 >>> Rational(1.5) 1 Rational can also accept strings that are valid literals for reals: >>> Rational("1.23") 123/100 >>> Rational('1e-2') 1/100 >>> Rational(".1") 1/10 Parsing needs for any other type of string for which a Rational is desired can be handled with the rational=True option in sympify() which produces rationals from strings like '.[3]' (=1/3) and '3/10' (=3/10). **Low-level** Access numerator and denominator as .p and .q: >>> r = Rational(3,4) >>> r 3/4 >>> r.p 3 >>> r.q 4 Note that p and q return integers (not sympy Integers) so some care is needed when using them in expressions: >>> r.p/r.q 0 """ is_real = True is_integer = False is_rational = True __slots__ = ['p', 'q'] is_Rational = True @cacheit def __new__(cls, p, q=None): if q is None: if isinstance(p, Rational): return p if isinstance(p, basestring): try: # we might have a Float neg_pow, digits, expt = decimal.Decimal(p).as_tuple() p = [1, -1][neg_pow] * int("".join(str(x) for x in digits)) if expt > 0: # TODO: this branch needs a test return Rational(p*Pow(10, expt), 1) return Rational(p, Pow(10, -expt)) except decimal.InvalidOperation: import re f = re.match('^([-+]?[0-9]+)/([0-9]+)$', p.replace(' ','')) if f: n, d = f.groups() return Rational(int(n), int(d)) raise ValueError('invalid literal: %s' % p) elif not isinstance(p, Basic): return Rational(S(p)) q = S.One if isinstance(q, Rational): p *= q.q q = q.p if isinstance(p, Rational): q *= p.q p = p.p p = int(p) q = int(q) if q == 0: if p == 0: if _errdict["divide"]: raise ValueError("Indeterminate 0/0") else: return S.NaN if p < 0: return S.NegativeInfinity return S.Infinity if q < 0: q = -q p = -p n = igcd(abs(p), q) if n > 1: p //= n q //= n if q == 1: return Integer(p) if p == 1 and q == 2: return S.Half obj = Expr.__new__(cls) obj.p = p obj.q = q #obj._args = (p, q) return obj
[docs] def limit_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. if max_denominator < 1: raise ValueError("max_denominator should be at least 1") if self.q <= max_denominator: return self p0, q0, p1, q1 = 0, 1, 1, 0 n, d = self.p, self.q while True: a = n//d q2 = q0+a*q1 if q2 > max_denominator: break p0, q0, p1, q1 = p1, q1, p0+a*p1, q2 n, d = d, n-a*d k = (max_denominator-q0)//q1 bound1 = Rational(p0+k*p1, q0+k*q1) bound2 = Rational(p1, q1) if abs(bound2 - self) <= abs(bound1-self): return bound2 else: return bound1
def __getnewargs__(self): return (self.p, self.q) def _hashable_content(self): return (self.p, self.q) def _eval_is_positive(self): return self.p > 0 def _eval_is_zero(self): return self.p == 0 def __neg__(self): return Rational(-self.p, self.q) @_sympifyit('other', NotImplemented) def __mul__(self, other): if (other is S.NaN) or (self is S.NaN): return S.NaN if isinstance(other, Float): return other * self if isinstance(other, Rational): return Rational(self.p * other.p, self.q * other.q) return Number.__mul__(self, other) @_sympifyit('other', NotImplemented) def __mod__(self, other): if isinstance(other, Rational): n = (self.p*other.q) // (other.p*self.q) return Rational(self.p*other.q - n*other.p*self.q, self.q*other.q) if isinstance(other, Float): return self.evalf() % other return Number.__mod__(self, other) @_sympifyit('other', NotImplemented) def __rmod__(self, other): if isinstance(other, Rational): return Rational.__mod__(other, self) if isinstance(other, Float): return other % self.evalf() return Number.__rmod__(self, other) # TODO reorder @_sympifyit('other', NotImplemented) def __add__(self, other): if (other is S.NaN) or (self is S.NaN): return S.NaN if isinstance(other, Float): return other + self if isinstance(other, Rational): if self.is_unbounded: if other.is_bounded: return self elif self==other: return self else: if other.is_unbounded: return other return Rational(self.p * other.q + self.q * other.p, self.q * other.q) return Number.__add__(self, other) def _eval_power(b, e): if (e is S.NaN): return S.NaN if isinstance(e, Number): if isinstance(e, Float): return b._eval_evalf(e._prec) ** e if e.is_negative: # (3/4)**-2 -> (4/3)**2 ne = -e if (ne is S.One): return Rational(b.q, b.p) if b < 0: if e.q != 1: return -(S.NegativeOne) ** ((e.p % e.q) / S(e.q)) * Rational(b.q, -b.p) ** ne else: return S.NegativeOne ** ne * Rational(b.q, -b.p) ** ne else: return Rational(b.q, b.p) ** ne if (e is S.Infinity): if b.p > b.q: # (3/2)**oo -> oo return S.Infinity if b.p < -b.q: # (-3/2)**oo -> oo + I*oo return S.Infinity + S.Infinity * S.ImaginaryUnit return S.Zero if isinstance(e, Integer): # (4/3)**2 -> 4**2 / 3**2 return Rational(b.p ** e.p, b.q ** e.p) if isinstance(e, Rational): if b.p != 1: # (4/3)**(5/6) -> 4**(5/6) * 3**(-5/6) return Integer(b.p) ** e * Integer(b.q) ** (-e) if b >= 0: return Integer(b.q)**Rational(e.p * (e.q-1), e.q) / ( Integer(b.q) ** Integer(e.p)) else: return (-1)**e * (-b)**e c, t = b.as_coeff_mul() if e.is_even and isinstance(c, Number) and c < 0: return (-c * Mul(*t)) ** e return def _as_mpf_val(self, prec): return mlib.from_rational(self.p, self.q, prec, rnd) def _mpmath_(self, prec, rnd): return mpmath.make_mpf(mlib.from_rational(self.p, self.q, prec, rnd)) def __abs__(self): return Rational(abs(self.p), self.q) def __int__(self): return int(float(self.p)/self.q) def __eq__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy != other --> not == if isinstance(other, NumberSymbol): if other.is_irrational: return False return other.__eq__(self) if isinstance(other, FunctionClass): #cos as opposed to cos(x) return False if other.is_comparable and not isinstance(other, Rational): other = other.evalf() if isinstance(other, Number): if isinstance(other, Float): return bool(mlib.mpf_eq(self._as_mpf_val(other._prec), other._mpf_)) return bool(self.p==other.p and self.q==other.q) return False # Rational != non-Number def __ne__(self, other): try: other = _sympify(other) except SympifyError: return True # sympy != other if isinstance(other, NumberSymbol): if other.is_irrational: return True return other.__ne__(self) if isinstance(other, FunctionClass): #cos as opposed to cos(x) return True if other.is_comparable and not isinstance(other, Rational): other = other.evalf() if isinstance(other, Number): if isinstance(other, Float): return bool(not mlib.mpf_eq(self._as_mpf_val(other._prec), other._mpf_)) return bool(self.p!=other.p or self.q!=other.q) return True # Rational != non-Number def __lt__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy > other --> not < if isinstance(other, NumberSymbol): return other.__ge__(self) if other.is_comparable and not isinstance(other, Rational): other = other.evalf() if isinstance(other, Number): if isinstance(other, Float): return bool(mlib.mpf_lt(self._as_mpf_val(other._prec), other._mpf_)) return bool(self.p * other.q < self.q * other.p) return Expr.__lt__(self, other) def __le__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy > other --> not <= if isinstance(other, NumberSymbol): return other.__gt__(self) if other.is_comparable and not isinstance(other, Rational): other = other.evalf() if isinstance(other, Number): if isinstance(other, Float): return bool(mlib.mpf_le(self._as_mpf_val(other._prec), other._mpf_)) return bool(self.p * other.q <= self.q * other.p) return Expr.__le__(self, other) def __hash__(self): return super(Rational, self).__hash__()
[docs] def factors(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. """ from sympy.ntheory import factorint f = factorint(self.p, limit=limit, use_trial=use_trial, use_rho=use_rho, use_pm1=use_pm1, verbose=verbose).copy() for p, e in factorint(self.q, limit=limit, use_trial=use_trial, use_rho=use_rho, use_pm1=use_pm1, verbose=verbose).items(): try: f[p] += -e except KeyError: f[p] = -e if len(f)>1 and 1 in f: del f[1] return f
[docs] def gcd(self, other): """Compute greatest common divisor of input arguments. """ if type(other) in (int, long): p = igcd(self.p, other) if self.is_Integer: return Integer(p) else: return Rational(p, self.q) else: other = _sympify(other) if other.is_Rational: p = igcd(self.p, other.p) if other.is_Integer: if self.is_Integer: return Integer(p) else: return Rational(p, self.q) else: if self.is_Integer: return Rational(p, other.q) else: return Rational(p, ilcm(self.q, other.q)) elif other.is_Number: return S.One else: raise TypeError("expected an integer or rational, got %s" % other)
[docs] def lcm(self, other): """Compute least common multiple of input arguments. """ if type(other) in (int, long): return Integer(ilcm(self.p, other)) else: other = _sympify(other) if other.is_Rational: p = ilcm(self.p, other.p) if self.is_Integer or other.is_Integer: return Integer(p) else: return Rational(p, igcd(self.q, other.q)) elif other.is_Number: return self*other else: raise TypeError("expected an integer or rational, got %s" % other)
[docs] def cofactors(self, other): """Compute GCD and cofactors of input arguments. """ other = _sympify(other) gcd = self.gcd(other) if gcd is S.One: return gcd, self, other else: return gcd, self/gcd, other/gcd
def as_numer_denom(self): return Integer(self.p), Integer(self.q) def _sage_(self): import sage.all as sage return sage.Integer(self.p)/sage.Integer(self.q) # int -> Integer
_intcache = {} # TODO move this tracing facility to sympy/core/trace.py ? def _intcache_printinfo(): ints = sorted(_intcache.keys()) nhit = _intcache_hits nmiss= _intcache_misses if nhit == 0 and nmiss == 0: print print 'Integer cache statistic was not collected' return miss_ratio = float(nmiss) / (nhit+nmiss) print print 'Integer cache statistic' print '-----------------------' print print '#items: %i' % len(ints) print print ' #hit #miss #total' print print '%5i %5i (%7.5f %%) %5i' % (nhit, nmiss, miss_ratio*100, nhit+nmiss) print print ints _intcache_hits = 0 _intcache_misses = 0 def int_trace(f): import os if os.getenv('SYMPY_TRACE_INT', 'no').lower() != 'yes': return f def Integer_tracer(cls, i): global _intcache_hits, _intcache_misses try: _intcache_hits += 1 return _intcache[i] except KeyError: _intcache_hits -= 1 _intcache_misses += 1 return f(cls, i) # also we want to hook our _intcache_printinfo into sys.atexit import atexit atexit.register(_intcache_printinfo) return Integer_tracer
[docs]class Integer(Rational): q = 1 is_integer = True is_Integer = True __slots__ = ['p'] def _as_mpf_val(self, prec): return mlib.from_int(self.p) def _mpmath_(self, prec, rnd): return mpmath.make_mpf(self._as_mpf_val(prec)) # TODO caching with decorator, but not to degrade performance @int_trace def __new__(cls, i): ival = int(i) try: return _intcache[ival] except KeyError: # We only work with well-behaved integer types. This converts, for # example, numpy.int32 instances. if ival == 0: obj = S.Zero elif ival == 1: obj = S.One elif ival == -1: obj = S.NegativeOne else: obj = Expr.__new__(cls) obj.p = ival _intcache[ival] = obj return obj def __getnewargs__(self): return (self.p,) # Arithmetic operations are here for efficiency def __int__(self): return self.p def __neg__(self): return Integer(-self.p) def __abs__(self): if self.p >= 0: return self else: return Integer(-self.p) def __divmod__(self, other): return divmod(self.p, other.p) # TODO make it decorator + bytecodehacks? def __add__(a, b): if isinstance(b, (int, long)): return Integer(a.p + b) elif isinstance(b, Integer): return Integer(a.p + b.p) return Rational.__add__(a, b) # a,b -not- b,a def __radd__(a, b): if isinstance(b, (int, long)): return Integer(b + a.p) elif isinstance(b, Integer): return Integer(b.p + a.p) return Rational.__add__(a, b) def __sub__(a, b): if isinstance(b, (int, long)): return Integer(a.p - b) elif isinstance(b, Integer): return Integer(a.p - b.p) return Rational.__sub__(a, b) def __rsub__(a, b): if isinstance(b, (int, long)): return Integer(b - a.p) elif isinstance(b, Integer): return Integer(b.p - a.p) return Rational.__rsub__(a, b) def __mul__(a, b): if isinstance(b, (int, long)): return Integer(a.p * b) elif isinstance(b, Integer): return Integer(a.p * b.p) return Rational.__mul__(a, b) def __rmul__(a, b): if isinstance(b, (int, long)): return Integer(b * a.p) elif isinstance(b, Integer): return Integer(b.p * a.p) return Rational.__mul__(a, b) def __mod__(a, b): if isinstance(b, (int, long)): return Integer(a.p % b) elif isinstance(b, Integer): return Integer(a.p % b.p) return Rational.__mod__(a, b) def __rmod__(a, b): if isinstance(b, (int, long)): return Integer(b % a.p) elif isinstance(b, Integer): return Integer(b.p % a.p) return Rational.__rmod__(a, b) def __eq__(a, b): if isinstance(b, (int, long)): return (a.p == b) elif isinstance(b, Integer): return (a.p == b.p) return Rational.__eq__(a, b) def __ne__(a, b): if isinstance(b, (int, long)): return (a.p != b) elif isinstance(b, Integer): return (a.p != b.p) return Rational.__ne__(a, b) def __gt__(a, b): if isinstance(b, (int, long)): return (a.p > b) elif isinstance(b, Integer): return (a.p > b.p) return Rational.__gt__(a, b) def __lt__(a, b): if isinstance(b, (int, long)): return (a.p < b) elif isinstance(b, Integer): return (a.p < b.p) return Rational.__lt__(a, b) def __ge__(a, b): if isinstance(b, (int, long)): return (a.p >= b) elif isinstance(b, Integer): return (a.p >= b.p) return Rational.__ge__(a, b) def __le__(a, b): if isinstance(b, (int, long)): return (a.p <= b) elif isinstance(b, Integer): return (a.p <= b.p) return Rational.__le__(a, b) def __hash__(self): return super(Integer, self).__hash__() def __index__(self): return self.p ######################################## def _eval_is_odd(self): return bool(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. """ from sympy import perfect_power if e is S.NaN: return S.NaN if b is S.One: return S.One if b is S.NegativeOne: return if e is S.Infinity: if b > S.One: return S.Infinity if b is S.NegativeOne: return S.NaN # cases for 0 and 1 are done in their respective classes return S.Infinity + S.ImaginaryUnit * S.Infinity if not isinstance(e, Number): # simplify when exp is even # (-2) ** k --> 2 ** k c, t = b.as_coeff_mul() if e.is_even and isinstance(c, Number) and c < 0: return (-c*Mul(*t))**e if not isinstance(e, Rational): return if e is S.Half and b < 0: # we extract I for this special case since everyone is doing so return S.ImaginaryUnit*Pow(-b, e) if e < 0: # invert base and change sign on exponent ne = -e if b < 0: if e.q != 1: return -(S.NegativeOne)**((e.p % e.q) / S(e.q)) * Rational(1, -b)**ne else: return (S.NegativeOne)**ne*Rational(1, -b)**ne else: return Rational(1, b)**ne # see if base is a perfect root, sqrt(4) --> 2 b_pos = int(abs(b)) x, xexact = integer_nthroot(b_pos, e.q) if xexact: # if it's a perfect root we've finished result = Integer(x ** abs(e.p)) if b < 0: result *= (-1)**e return result # 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 power p = perfect_power(b_pos) if p is not False: dict = {p[0]: p[1]} else: dict = Integer(b_pos).factors(limit=2**15) # now process the dict of factors if b.is_negative: dict[-1] = 1 out_int = 1 # integer part out_rad = 1 # extracted radicals sqr_int = 1 sqr_gcd = 0 sqr_dict = {} for prime, exponent in dict.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) if div_e > 0: out_int *= prime**div_e if div_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) if g != 1: out_rad *= Pow(prime, Rational(div_m//g, e.q//g)) else: sqr_dict[prime] = div_m # identify gcd of remaining powers for p, ex in sqr_dict.iteritems(): if sqr_gcd == 0: sqr_gcd = ex else: sqr_gcd = igcd(sqr_gcd, ex) if sqr_gcd == 1: break for k, v in sqr_dict.iteritems(): sqr_int *= k**(v//sqr_gcd) if sqr_int == b and out_int == 1 and out_rad == 1: result = None else: result = out_int*out_rad*Pow(sqr_int, Rational(sqr_gcd, e.q)) return result def _eval_is_prime(self): if self.p < 0: return False def as_numer_denom(self): return self, S.One def __floordiv__(self, other): return Integer(self.p // Integer(other).p) def __rfloordiv__(self, other): return Integer(Integer(other).p // self.p)
[docs] def factorial(a): """Compute factorial of `a`. """ from sympy.functions.combinatorial.factorials import factorial return Integer(factorial(int(a)))
[docs] def isqrt(a): """Compute integer square root of `a`. """ return Integer(mlib.isqrt(int(a)))
[docs] def half_gcdex(a, b): """Half Extended Euclidean Algorithm. """ s, _, h = a.gcdex(b) return s, h
[docs] def gcdex(a, b): """Extended Euclidean Algorithm. """ if isinstance(b, (int, long)): return tuple(map(Integer, igcdex(int(a), b))) else: b = _sympify(b) if b.is_Integer: return tuple(map(Integer, igcdex(int(a), int(b)))) else: raise ValueError("expected an integer, got %s" % b)
[docs] def invert(a, b): """Invert `a` modulo `b`, if possible. """ if isinstance(b, (int, long)): a = int(a) else: b = _sympify(b) if b.is_Integer: a, b = int(a), int(b) else: raise ValueError("expected an integer, got %s" % b) s, _, h = igcdex(a, b) if h == 1: return Integer(s % b) else: raise ZeroDivisionError("zero divisor") # Add sympify converters
converter[int] = converter[long] = Integer class RationalConstant(Rational): """ Abstract base class for rationals with specific behaviors Derived classes must define class attributes p and q and should probably all be singletons. """ __slots__ = [] def __new__(cls): return AtomicExpr.__new__(cls) class IntegerConstant(Integer): __slots__ = [] def __new__(cls): return AtomicExpr.__new__(cls) class Zero(IntegerConstant): __metaclass__ = Singleton p = 0 q = 1 is_positive = False is_negative = False is_finite = False is_zero = True is_prime = False is_composite = False __slots__ = [] @staticmethod def __abs__(): return S.Zero @staticmethod def __neg__(): return S.Zero def _eval_power(b, e): if e.is_negative: return S.Infinity if e.is_positive: return b d = e.evalf() if isinstance(d, Number): if d.is_negative: return S.Infinity return b coeff, terms = e.as_coeff_mul() if coeff.is_negative: return S.Infinity ** Mul(*terms) if coeff is not S.One: return b ** Mul(*terms) def _eval_order(self, *symbols): # Order(0,x) -> 0 return self def __nonzero__(self): return False class One(IntegerConstant): __metaclass__ = Singleton p = 1 q = 1 is_prime = True __slots__ = [] def _eval_evalf(self, prec): return self @staticmethod def __abs__(): return S.One @staticmethod def __neg__(): return S.NegativeOne def _eval_order(self, *symbols): return @staticmethod def factors(): return {1: 1} class NegativeOne(IntegerConstant): __metaclass__ = Singleton p = -1 q = 1 __slots__ = [] def _eval_evalf(self, prec): return self @staticmethod def __abs__(): return S.One @staticmethod def __neg__(): return S.One def _eval_power(b, e): if e.is_odd: return S.NegativeOne if e.is_even: return S.One if isinstance(e, Number): if isinstance(e, Float): return Float(-1.0) ** e if e is S.NaN: return S.NaN if e is S.Infinity or e is S.NegativeInfinity: return S.NaN if e is S.Half: return S.ImaginaryUnit if isinstance(e, Rational): if e.q == 2: return S.ImaginaryUnit ** Integer(e.p) q = Float(e).floor() if q: q = Integer(q) return b ** q * b ** (e - q) return class Half(RationalConstant): __metaclass__ = Singleton p = 1 q = 2 __slots__ = [] @staticmethod def __abs__(): return S.Half class Infinity(RationalConstant): __metaclass__ = Singleton p = 1 q = 0 __slots__ = [] is_commutative = True is_positive = True is_bounded = False is_finite = False is_infinitesimal = False is_integer = None is_rational = None is_odd = None @staticmethod def __abs__(): return S.Infinity @staticmethod def __neg__(): return S.NegativeInfinity def _eval_power(b, e): """ e is symbolic object but not equal to 0, 1 oo ** nan -> nan oo ** (-p) -> 0, p is number, oo """ if e.is_positive: return S.Infinity if e.is_negative: return S.Zero if isinstance(e, Number): if e is S.NaN: return S.NaN d = e.evalf() if isinstance(d, Number): return b ** d return def _as_mpf_val(self, prec): return mlib.finf def _sage_(self): import sage.all as sage return sage.oo def __gt__(a, b): if b is S.Infinity: return False return True def __lt__(a, b): return False def __ge__(a, b): return True def __le__(a, b): if b is S.Infinity: return True return False def __mod__(self, other): return S.NaN __rmod__ = __mod__ oo = S.Infinity class NegativeInfinity(RationalConstant): __metaclass__ = Singleton p = -1 q = 0 __slots__ = [] is_commutative = True is_real = True is_positive = False is_bounded = False is_finite = False is_infinitesimal = False is_integer = None is_rational = None @staticmethod def __abs__(): return S.Infinity @staticmethod def __neg__(): return S.Infinity def _eval_power(b, e): """ e is symbolic object but not equal to 0, 1 (-oo) ** nan -> nan (-oo) ** oo -> nan (-oo) ** (-oo) -> nan (-oo) ** e -> oo, e is positive even integer (-oo) ** o -> -oo, o is positive odd integer """ if isinstance(e, Number): if (e is S.NaN) or (e is S.Infinity) or (e is S.NegativeInfinity): return S.NaN if isinstance(e, Integer): if e.is_positive: if e.is_odd: return S.NegativeInfinity return S.Infinity return S.NegativeOne**e * S.Infinity ** e return def _as_mpf_val(self, prec): return mlib.fninf def _sage_(self): import sage.all as sage return -(sage.oo) def __gt__(a, b): return False def __lt__(a, b): if b is S.NegativeInfinity: return False return True def __ge__(a, b): if b is S.NegativeInfinity: return True return False def __le__(a, b): return True class NaN(RationalConstant): __metaclass__ = Singleton p = 0 q = 0 is_commutative = True is_real = None is_rational = None is_integer = None is_comparable = False is_finite = None is_bounded = None #is_unbounded = False is_zero = None is_prime = None is_positive = None __slots__ = [] def _as_mpf_val(self, prec): return mlib.fnan def _eval_power(b, e): if e is S.Zero: return S.One return b def _sage_(self): import sage.all as sage return sage.NaN nan = S.NaN class ComplexInfinity(AtomicExpr): __metaclass__ = Singleton is_commutative = True is_comparable = None is_bounded = False is_real = None is_number = True __slots__ = [] def __new__(cls): return AtomicExpr.__new__(cls) @staticmethod def __abs__(): return S.Infinity @staticmethod def __neg__(): return S.ComplexInfinity def _eval_power(b, e): if e is S.ComplexInfinity: return S.NaN if isinstance(e, Number): if e is S.Zero: return S.NaN else: if e.is_positive: return S.ComplexInfinity else: return S.Zero zoo = S.ComplexInfinity
[docs]class NumberSymbol(AtomicExpr): __metaclass__ = Singleton is_commutative = True is_comparable = True is_bounded = True is_finite = True is_number = True __slots__ = [] is_NumberSymbol = True def __new__(cls): return AtomicExpr.__new__(cls)
[docs] def approximation(self, number_cls): """ Return an interval with number_cls endpoints that contains the value of NumberSymbol. If not implemented, then return None. """
def _eval_evalf(self, prec): return Float._new(self._as_mpf_val(prec), prec) def __eq__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy != other --> not == if self is other: return True if isinstance(other, Number) and self.is_irrational: return False return False # NumberSymbol != non-(Number|self) def __ne__(self, other): try: other = _sympify(other) except SympifyError: return True # sympy != other if self is other: return False if isinstance(other, Number) and self.is_irrational: return True return True # NumberSymbol != non(Number|self) def __lt__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy > other --> not < if self is other: return False if isinstance(other, Number): approx = self.approximation_interval(other.__class__) if approx is not None: l,u = approx if other < l: return False if other > u: return True return self.evalf()<other if other.is_comparable: other = other.evalf() return self.evalf()<other return Expr.__lt__(self, other) def __le__(self, other): try: other = _sympify(other) except SympifyError: return False # sympy > other --> not <= if self is other: return True if other.is_comparable: other = other.evalf() if isinstance(other, Number): return self.evalf()<=other return Expr.__le__(self, other) def __gt__(self, other): return (-self) < (-other) def __ge__(self, other): return (-self) <= (-other) def __int__(self): return int(self.evalf(0)) def __hash__(self): return super(NumberSymbol, self).__hash__()
class Exp1(NumberSymbol): __metaclass__ = Singleton is_real = True is_positive = True is_negative = False # XXX Forces is_negative/is_nonnegative is_irrational = True __slots__ = [] @staticmethod def __abs__(): return S.Exp1 def _as_mpf_val(self, prec): return mpf_e(prec) def approximation_interval(self, number_cls): if issubclass(number_cls,Integer): return (Integer(2),Integer(3)) elif issubclass(number_cls,Rational): pass def _eval_power(self, exp): return C.exp(exp) def _sage_(self): import sage.all as sage return sage.e E = S.Exp1 class Pi(NumberSymbol): __metaclass__ = Singleton is_real = True is_positive = True is_negative = False is_irrational = True __slots__ = [] @staticmethod def __abs__(): return S.Pi def _as_mpf_val(self, prec): return mpf_pi(prec) def approximation_interval(self, number_cls): if issubclass(number_cls, Integer): return (Integer(3), Integer(4)) elif issubclass(number_cls, Rational): return (Rational(223,71), Rational(22,7)) def _sage_(self): import sage.all as sage return sage.pi pi = S.Pi class GoldenRatio(NumberSymbol): __metaclass__ = Singleton is_real = True is_positive = True is_negative = False is_irrational = True __slots__ = [] def _as_mpf_val(self, prec): return mlib.from_man_exp(phi_fixed(prec+10), -prec-10) def _eval_expand_func(self, deep=True, **hints): return S.Half + S.Half*S.Sqrt(5) def approximation_interval(self, number_cls): if issubclass(number_cls, Integer): return (S.One, Rational(2)) elif issubclass(number_cls, Rational): pass def _sage_(self): import sage.all as sage return sage.golden_ratio class EulerGamma(NumberSymbol): __metaclass__ = Singleton is_real = True is_positive = True is_negative = False is_irrational = None __slots__ = [] def _as_mpf_val(self, prec): return mlib.from_man_exp(mlib.libhyper.euler_fixed( prec+10), -prec-10) def approximation_interval(self, number_cls): if issubclass(number_cls, Integer): return (S.Zero, S.One) elif issubclass(number_cls, Rational): return (S.Half, Rational(3, 5)) def _sage_(self): import sage.all as sage return sage.euler_gamma class Catalan(NumberSymbol): __metaclass__ = Singleton is_real = True is_positive = True is_negative = False is_irrational = None __slots__ = [] def _as_mpf_val(self, prec): return mlib.from_man_exp(mlib.catalan_fixed(prec+10), -prec-10) def approximation_interval(self, number_cls): if issubclass(number_cls, Integer): return (S.Zero, S.One) elif issubclass(number_cls, Rational): return (Rational(9, 10), S.One) def _sage_(self): import sage.all as sage return sage.catalan class ImaginaryUnit(AtomicExpr): __metaclass__ = Singleton is_commutative = True is_imaginary = True is_bounded = True is_finite = True is_number = True __slots__ = [] @staticmethod def __abs__(): return S.One def _eval_evalf(self, prec): return self def _eval_conjugate(self): return -S.ImaginaryUnit def _eval_power(b, e): """ b is I = sqrt(-1) e is symbolic object but not equal to 0, 1 I ** r -> (-1)**(r/2) -> exp(r/2 * Pi * I) -> sin(Pi*r/2) + cos(Pi*r/2) * I, r is decimal I ** 0 mod 4 -> 1 I ** 1 mod 4 -> I I ** 2 mod 4 -> -1 I ** 3 mod 4 -> -I """ if isinstance(e, Number): if isinstance(e, Integer): ei = e.p % 4 if ei == 0: return S.One if ei == 1: return S.ImaginaryUnit if ei == 2: return -S.One return -S.ImaginaryUnit return (S.NegativeOne) ** (e * S.Half) return def as_base_exp(self): return S.NegativeOne, S.Half def _sage_(self): import sage.all as sage return sage.I I = S.ImaginaryUnit try: # fractions is only available for python 2.6+ import fractions def sympify_fractions(f): return Rational(f.numerator, f.denominator) converter[fractions.Fraction] = sympify_fractions except ImportError: pass try: import gmpy def sympify_mpz(x): return Integer(long(x)) def sympify_mpq(x): return Rational(long(x.numer()), long(x.denom())) converter[type(gmpy.mpz(1))] = sympify_mpz converter[type(gmpy.mpq(1, 2))] = sympify_mpq except ImportError: pass def sympify_mpmath(x): return Expr._from_mpmath(x, x.context.prec) converter[mpnumeric] = sympify_mpmath def sympify_complex(a): real, imag = map(sympify, (a.real, a.imag)) return real + S.ImaginaryUnit * imag converter[complex] = sympify_complex _intcache[0] = S.Zero _intcache[1] = S.One _intcache[-1]= S.NegativeOne from function import FunctionClass from power import Pow, integer_nthroot from mul import Mul Mul.identity = One() from add import Add Add.identity = Zero()