sympy/core/numbers.py

   1 import math
   2 import sympy.mpmath as mpmath
   3 import sympy.mpmath.lib as mlib
   4 import sympy.mpmath.libmpc as mlibc
   5 import decimal
   6
   7 rnd = mlib.round_nearest
   8
   9 from basic import Basic, Atom, S, C, SingletonMeta, Memoizer, MemoizerArg
  10 from sympify import _sympify, SympifyError, _sympifyit
  11 from power import integer_nthroot
  12
  13 # from mul import Mul   /cyclic/
  14 # from power import Pow /cyclic/
  15 # from function import FunctionClass    /cyclic/
  16
  17 # (a,b) -> gcd(a,b)
  18 _gcdcache = {}
  19
  20 # TODO caching with decorator, but not to degrade performance
  21 def igcd(a, b):
  22     """Computes integer greates common divisor of two numbers.
  23
  24        The algorithm is based on the well known Euclid's algorithm. To
  25        improve speed igcd() has its own caching mechanizm implemented.
  26     """
  27     try:
  28         return _gcdcache[(a,b)]
  29     except KeyError:
  30         if a and b:
  31             if b < 0:
  32                 b = -b
  33
  34             while b:
  35                 a, b = b, a % b
  36         else:
  37             a = abs(a or b)
  38
  39         _gcdcache[(a,b)] = a
  40         return a
  41
  42 def ilcm(a, b):
  43     """Computes integer least common multiple of two numbers. """
  44     if a == 0 and b == 0:
  45         return 0
  46     else:
  47         return a * b / igcd(a, b)
  48
  49 def igcdex(a, b):
  50     """Returns x, y, g such that g = x*a + y*b = gcd(a, b).
  51
  52        >>> igcdex(2, 3)
  53        (-1, 1, 1)
  54        >>> igcdex(10, 12)
  55        (-1, 1, 2)
  56
  57        >>> x, y, g = igcdex(100, 2004)
  58        >>> x, y, g
  59        (-20, 1, 4)
  60        >>> x*100 + y*2004
  61        4
  62
  63     """
  64     if (not a) and (not b):
  65         return (0, 1, 0)
  66
  67     if not a:
  68         return (0, b/abs(b), abs(b))
  69     if not b:
  70         return (a/abs(a), 0, abs(a))
  71
  72     if a < 0:
  73         a, x_sign = -a, -1
  74     else:
  75         x_sign = 1
  76
  77     if b < 0:
  78         b, y_sign = -b, -1
  79     else:
  80         y_sign = 1
  81
  82     x, y, r, s = 1, 0, 0, 1
  83
  84     while b:
  85         (c, q) = (a % b, a / b)
  86         (a, b, r, s, x, y) = (b, c, x-q*r, y-q*s, r, s)
  87
  88     return (x*x_sign, y*y_sign, a)
  89
  90 @Memoizer((int, long), return_value_converter = lambda d: d.copy())
  91 def factor_trial_division(n):
  92     """
  93     Factor any integer into a product of primes, 0, 1, and -1.
  94     Returns a dictionary {<prime: exponent>}.
  95     """
  96     if not n:
  97         return {0:1}
  98     factors = {}
  99     if n < 0:
 100         factors[-1] = 1
 101         n = -n
 102     if n==1:
 103         factors[1] = 1
 104         return factors
 105     d = 2
 106     while n % d == 0:
 107         try:
 108             factors[d] += 1
 109         except KeyError:
 110             factors[d] = 1
 111         n //= d
 112     d = 3
 113     while n > 1 and d*d <= n:
 114         if n % d:
 115             d += 2
 116         else:
 117             try:
 118                 factors[d] += 1
 119             except KeyError:
 120                 factors[d] = 1
 121             n //= d
 122     if n>1:
 123         try:
 124             factors[n] += 1
 125         except KeyError:
 126             factors[n] = 1
 127     return factors
 128
 129
 130 class Number(Atom):
 131     """
 132     Represents any kind of number in sympy.
 133
 134     Floating point numbers are represented by the Real class.
 135     Integer numbers (of any size), together with rational numbers (again, there
 136     is no limit on their size) are represented by the Rational class.
 137
 138     If you want to represent for example 1+sqrt(2), then you need to do:
 139
 140     Rational(1) + sqrt(Rational(2))
 141     """
 142     is_commutative = True
 143     is_comparable = True
 144     is_bounded = True
 145     is_finite = True
 146
 147     __slots__ = []
 148
 149     # Used to make max(x._prec, y._prec) return x._prec when only x is a float
 150     _prec = -1
 151
 152     is_Number = True
 153
 154     def __new__(cls, *obj):
 155         if len(obj)==1: obj=obj[0]
 156         if isinstance(obj, (int, long)):
 157             return Integer(obj)
 158         if isinstance(obj,tuple) and len(obj)==2:
 159             return Rational(*obj)
 160         if isinstance(obj, (str,float,mpmath.mpf,decimal.Decimal)):
 161             return Real(obj)
 162         if isinstance(obj, Number):
 163             return obj
 164         raise TypeError("expected str|int|long|float|Decimal|Number object but got %r" % (obj))
 165
 166     def _as_mpf_val(self, prec):
 167         """Evaluate to mpf tuple accurate to at least prec bits"""
 168         raise NotImplementedError('%s needs ._as_mpf_val() method' % \
 169             (self.__class__.__name__))
 170
 171     def _eval_evalf(self, prec):
 172         return Real._new(self._as_mpf_val(prec), prec)
 173
 174     def _as_mpf_op(self, prec):
 175         prec = max(prec, self._prec)
 176         return self._as_mpf_val(prec), prec
 177
 178     def __float__(self):
 179         return mlib.to_float(self._as_mpf_val(53))
 180
 181     def _eval_derivative(self, s):
 182         return S.Zero
 183
 184     def _eval_conjugate(self):
 185         return self
 186
 187     def _eval_order(self, *symbols):
 188         # Order(5, x, y) -> Order(1,x,y)
 189         return C.Order(S.One, *symbols)
 190
 191     def __eq__(self, other):
 192         raise NotImplementedError('%s needs .__eq__() method' % (self.__class__.__name__))
 193     def __ne__(self, other):
 194         raise NotImplementedError('%s needs .__ne__() method' % (self.__class__.__name__))
 195     def __lt__(self, other):
 196         raise NotImplementedError('%s needs .__lt__() method' % (self.__class__.__name__))
 197     def __le__(self, other):
 198         raise NotImplementedError('%s needs .__le__() method' % (self.__class__.__name__))
 199
 200     def __gt__(self, other):
 201         return _sympify(other).__lt__(self)
 202     def __ge__(self, other):
 203         return _sympify(other).__le__(self)
 204
 205     def as_coeff_terms(self, x=None):
 206         # a -> c * t
 207         return self, tuple()
 208
 209
 210 class Real(Number):
 211     """
 212     Represents a floating point number. It is capable of representing
 213     arbitrary-precision floating-point numbers
 214
 215     Usage:
 216     ======
 217         Real(3.5)   .... 3.5 (the 3.5 was converted from a python float)
 218         Real("3.0000000000000005")
 219
 220     Notes:
 221     ======
 222         - Real(x) with x being a Python int/long will return Integer(x)
 223     """
 224     is_real = True
 225     is_irrational = False
 226     is_integer = False
 227
 228     __slots__ = ['_mpf_', '_prec']
 229
 230     # mpz can't be pickled
 231     def __getstate__(self):
 232         d = Basic.__getstate__(self).copy()
 233         del d["_mpf_"]
 234         return mlib.to_pickable(self._mpf_), d
 235
 236     def __setstate__(self, state):
 237         _mpf_, d = state
 238         _mpf_ = mlib.from_pickable(_mpf_)
 239         self._mpf_ = _mpf_
 240         Basic.__setstate__(self, d)
 241
 242     is_Real = True
 243
 244     def floor(self):
 245         return C.Integer(int(mlib.to_int(mlib.ffloor(self._mpf_, self._prec))))
 246
 247     def ceiling(self):
 248         return C.Integer(int(mlib.to_int(mlib.fceil(self._mpf_, self._prec))))
 249
 250     @property
 251     def num(self):
 252         return mpmath.mpf(self._mpf_)
 253
 254     def _as_mpf_val(self, prec):
 255         return self._mpf_
 256
 257     def _as_mpf_op(self, prec):
 258         return self._mpf_, max(prec, self._prec)
 259
 260     def __new__(cls, num, prec=15):
 261         prec = mlib.dps_to_prec(prec)
 262         if isinstance(num, (int, long)):
 263             return Integer(num)
 264         if isinstance(num, (str, decimal.Decimal)):
 265             _mpf_ = mlib.from_str(str(num), prec, rnd)
 266         elif isinstance(num, tuple) and len(num) == 4:
 267             _mpf_ = num
 268         else:
 269             _mpf_ = mpmath.mpf(num)._mpf_
 270         if not num:
 271             return C.Zero()
 272         obj = Basic.__new__(cls)
 273         obj._mpf_ = _mpf_
 274         obj._prec = prec
 275         return obj
 276
 277     @classmethod
 278     def _new(cls, _mpf_, _prec):
 279         if _mpf_ == mlib.fzero:
 280             return S.Zero
 281         obj = Basic.__new__(cls)
 282         obj._mpf_ = _mpf_
 283         obj._prec = _prec
 284         return obj
 285
 286     def _hashable_content(self):
 287         return (self._mpf_, self._prec)
 288
 289     def _eval_is_positive(self):
 290         return self.num > 0
 291
 292     def _eval_is_negative(self):
 293         return self.num < 0
 294
 295     def __neg__(self):
 296         return Real._new(mlib.fneg(self._mpf_), self._prec)
 297
 298     def __mul__(self, other):
 299         try:
 300             other = _sympify(other)
 301         except SympifyError:
 302             return NotImplemented
 303         if isinstance(other, Number):
 304             rhs, prec = other._as_mpf_op(self._prec)
 305             return Real._new(mlib.fmul(self._mpf_, rhs, prec, rnd), prec)
 306         return Number.__mul__(self, other)
 307
 308     def __add__(self, other):
 309         try:
 310             other = _sympify(other)
 311         except SympifyError:
 312             return NotImplemented
 313         if (other is S.NaN) or (self is NaN):
 314             return S.NaN
 315         if isinstance(other, Number):
 316             rhs, prec = other._as_mpf_op(self._prec)
 317             return Real._new(mlib.fadd(self._mpf_, rhs, prec, rnd), prec)
 318         return Number.__add__(self, other)
 319
 320     def _eval_power(b, e):
 321         """
 322         b is Real but not equal to rationals, integers, 0.5, oo, -oo, nan
 323         e is symbolic object but not equal to 0, 1
 324
 325         (-p) ** r -> exp(r * log(-p)) -> exp(r * (log(p) + I*Pi)) ->
 326                   -> p ** r * (sin(Pi*r) + cos(Pi*r) * I)
 327         """
 328         if isinstance(e, Number):
 329             if isinstance(e, Integer):
 330                 prec = b._prec
 331                 return Real._new(mlib.fpowi(b._mpf_, e.p, prec, rnd), prec)
 332             e, prec = e._as_mpf_op(b._prec)
 333             b = b._mpf_
 334             try:
 335                 y = mlib.fpow(b, e, prec, rnd)
 336                 return Real._new(y, prec)
 337             except mlib.ComplexResult:
 338                 re, im = mlibc.mpc_pow((b, mlib.fzero), (e, mlib.fzero), prec, rnd)
 339                 return Real._new(re, prec) + Real._new(im, prec) * S.ImaginaryUnit
 340
 341     def __abs__(self):
 342         return Real._new(mlib.fabs(self._mpf_), self._prec)
 343
 344     def __int__(self):
 345         return int(mlib.to_int(self._mpf_))
 346
 347     def __eq__(self, other):
 348         try:
 349             other = _sympify(other)
 350         except SympifyError:
 351             return False    # sympy != other  -->  not ==
 352         if isinstance(other, NumberSymbol):
 353             if other.is_irrational: return False
 354             return other.__eq__(self)
 355         if other.is_comparable: other = other.evalf()
 356         if isinstance(other, Number):
 357             return bool(mlib.feq(self._mpf_, other._as_mpf_val(self._prec)))
 358         return False    # Real != non-Number
 359
 360     def __ne__(self, other):
 361         try:
 362             other = _sympify(other)
 363         except SympifyError:
 364             return True     # sympy != other
 365         if isinstance(other, NumberSymbol):
 366             if other.is_irrational: return True
 367             return other.__ne__(self)
 368         if other.is_comparable: other = other.evalf()
 369         if isinstance(other, Number):
 370             return bool(not mlib.feq(self._mpf_, other._as_mpf_val(self._prec)))
 371         return True     # Real != non-Number
 372
 373     def __lt__(self, other):
 374         try:
 375             other = _sympify(other)
 376         except SympifyError:
 377             return False    # sympy > other
 378         if isinstance(other, NumberSymbol):
 379             return other.__ge__(self)
 380         if other.is_comparable: other = other.evalf()
 381         if isinstance(other, Number):
 382             return bool(mlib.flt(self._mpf_, other._as_mpf_val(self._prec)))
 383         return Basic.__lt__(self, other)
 384
 385     def __le__(self, other):
 386         try:
 387             other = _sympify(other)
 388         except SympifyError:
 389             return False    # sympy > other  -->  ! <=
 390         if isinstance(other, NumberSymbol):
 391             return other.__gt__(self)
 392         if other.is_comparable: other = other.evalf()
 393         if isinstance(other, Number):
 394             return bool(mlib.fle(self._mpf_, other._as_mpf_val(self._prec)))
 395         return Basic.__le__(self, other)
 396
 397     def epsilon_eq(self, other, epsilon="10e-16"):
 398         return abs(self - other) < Real(epsilon)
 399
 400 # this is here to work nicely in Sage
 401 RealNumber = Real
 402
 403
 404 def _parse_rational(s):
 405     """Parse rational number from string representation"""
 406     # Simple fraction
 407     if "/" in s:
 408         p, q = s.split("/")
 409         return int(p), int(q)
 410     # Recurring decimal
 411     elif "[" in s:
 412         sign = 1
 413         if s[0] == "-":
 414             sign = -1
 415             s = s[1:]
 416         s, periodic = s.split("[")
 417         periodic = periodic.rstrip("]")
 418         offset = len(s) - s.index(".") - 1
 419         n1 = int(periodic)
 420         n2 = int("9" * len(periodic))
 421         r = Rational(*_parse_rational(s)) + Rational(n1, n2*10**offset)
 422         return sign*r.p, r.q
 423     # Ordinary decimal string. Use the Decimal class's built-in parser
 424     else:
 425         sign, digits, expt = decimal.Decimal(s).as_tuple()
 426         p = (1, -1)[sign] * int("".join(str(x) for x in digits))
 427         if expt >= 0:
 428             return p*(10**expt), 1
 429         else:
 430             return p, 10**-expt
 431
 432
 433 class Rational(Number):
 434     """Represents integers and rational numbers (p/q) of any size.
 435
 436     Examples
 437     ========
 438         >>> Rational(3)
 439         3
 440         >>> Rational(1,2)
 441         1/2
 442
 443     You can create a rational from a string:
 444         >>> Rational("3/5")
 445         3/5
 446         >>> Rational("1.23")
 447         123/100
 448
 449     Use square brackets to indicate a recurring decimal:
 450         >>> Rational("0.[333]")
 451         1/3
 452         >>> Rational("1.2[05]")
 453         1193/990
 454         >>> float(Rational(1193,990))
 455         1.20505050505
 456
 457
 458     Low-level
 459     ---------
 460
 461     Access nominator and denominator as .p and .q:
 462         >>> r = Rational(3,4)
 463         >>> r
 464         3/4
 465         >>> r.p
 466         3
 467         >>> r.q
 468         4
 469     """
 470     is_real = True
 471     is_integer = False
 472     is_rational = True
 473
 474     __slots__ = ['p', 'q']
 475
 476     is_Rational = True
 477
 478     @Memoizer(type, (int, long, str), MemoizerArg((int, long, type(None)), name="q"))
 479     def __new__(cls, p, q = None):
 480         if q is None:
 481             if isinstance(p, str):
 482                 p, q = _parse_rational(p)
 483             else:
 484                 return Integer(p)
 485         if q==0:
 486             if p==0: return S.NaN
 487             if p<0: return S.NegativeInfinity
 488             return S.Infinity
 489         if q<0:
 490             q = -q
 491             p = -p
 492         n = igcd(abs(p), q)
 493         if n>1:
 494             p /= n
 495             q /= n
 496         if q==1: return Integer(p)
 497         if p==1 and q==2: return S.Half
 498         obj = Basic.__new__(cls)
 499         obj.p = p
 500         obj.q = q
 501         #obj._args = (p, q)
 502         return obj
 503
 504     def _hashable_content(self):
 505         return (self.p, self.q)
 506
 507     def _eval_is_positive(self):
 508         return self.p > 0
 509
 510     def _eval_is_zero(self):
 511         return self.p == 0
 512
 513     def __neg__(self): return Rational(-self.p, self.q)
 514
 515     @_sympifyit('other', NotImplemented)
 516     def __mul__(self, other):
 517         if (other is S.NaN) or (self is S.NaN):
 518             return S.NaN
 519         if isinstance(other, Real):
 520             return other * self
 521         if isinstance(other, Rational):
 522             return Rational(self.p * other.p, self.q * other.q)
 523         return Number.__mul__(self, other)
 524
 525     # TODO reorder
 526     @_sympifyit('other', NotImplemented)
 527     def __add__(self, other):
 528         if (other is S.NaN) or (self is S.NaN):
 529             return S.NaN
 530         if isinstance(other, Real):
 531             return other + self
 532         if isinstance(other, Rational):
 533             if self.is_unbounded:
 534                 if other.is_bounded:
 535                     return self
 536                 elif self==other:
 537                     return self
 538             else:
 539                 if other.is_unbounded:
 540                     return other
 541             return Rational(self.p * other.q + self.q * other.p, self.q * other.q)
 542         return Number.__add__(self, other)
 543
 544     def _eval_power(b, e):
 545         if isinstance(e, Number):
 546             if (e is S.NaN): return S.NaN
 547             if isinstance(e, Real):
 548                 return b._eval_evalf(e._prec) ** e
 549             if e.is_negative:
 550                 # (3/4)**-2 -> (4/3)**2
 551                 ne = -e
 552                 if (ne is S.One):
 553                     return Rational(b.q, b.p)
 554                 return Rational(b.q, b.p) ** ne
 555             if (e is S.Infinity):
 556                 if b.p > b.q:
 557                     # (3/2)**oo -> oo
 558                     return S.Infinity
 559                 if b.p < -b.q:
 560                     # (-3/2)**oo -> oo + I*oo
 561                     return S.Infinity + S.Infinity * S.ImaginaryUnit
 562                 return S.Zero
 563             if isinstance(e, Integer):
 564                 # (4/3)**2 -> 4**2 / 3**2
 565                 return Rational(b.p ** e.p, b.q ** e.p)
 566             if isinstance(e, Rational):
 567                 if b.p != 1:
 568                     # (4/3)**(5/6) -> 4**(5/6) * 3**(-5/6)
 569                     return Integer(b.p) ** e * Integer(b.q) ** (-e)
 570                 if b >= 0:
 571                     return Integer(b.q)**Rational(e.p * (e.q-1), e.q) / ( Integer(b.q) ** Integer(e.p))
 572                 else:
 573                     return (-1)**e * (-b)**e
 574
 575         c,t = b.as_coeff_terms()
 576         if e.is_even and isinstance(c, Number) and c < 0:
 577             return (-c * Mul(*t)) ** e
 578
 579         return
 580
 581     def _as_mpf_val(self, prec):
 582         return mlib.from_rational(self.p, self.q, prec, rnd)
 583
 584     def __abs__(self):
 585         return Rational(abs(self.p), self.q)
 586
 587     def __int__(self):
 588         return int(self.p//self.q)
 589
 590     def __eq__(self, other):
 591         try:
 592             other = _sympify(other)
 593         except SympifyError:
 594             return False    # sympy != other  -->  not ==
 595         if isinstance(other, NumberSymbol):
 596             if other.is_irrational: return False
 597             return other.__eq__(self)
 598         if isinstance(self, Number) and isinstance(other, FunctionClass):
 599             return False
 600         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 601         if isinstance(other, Number):
 602             if isinstance(other, Real):
 603                 return bool(mlib.feq(self._as_mpf_val(other._prec), other._mpf_))
 604             return bool(self.p==other.p and self.q==other.q)
 605
 606         return False    # Rational != non-Number
 607
 608     def __ne__(self, other):
 609         try:
 610             other = _sympify(other)
 611         except SympifyError:
 612             return True     # sympy != other
 613         if isinstance(other, NumberSymbol):
 614             if other.is_irrational: return True
 615             return other.__ne__(self)
 616         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 617         if isinstance(other, Number):
 618             if isinstance(other, Real):
 619                 return bool(not mlib.feq(self._as_mpf_val(other._prec), other._mpf_))
 620             return bool(self.p!=other.p or self.q!=other.q)
 621
 622         return True     # Rational != non-Number
 623
 624     def __lt__(self, other):
 625         try:
 626             other = _sympify(other)
 627         except SympifyError:
 628             return False    # sympy > other  --> not <
 629         if isinstance(other, NumberSymbol):
 630             return other.__ge__(self)
 631         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 632         if isinstance(other, Number):
 633             if isinstance(other, Real):
 634                 return bool(mlib.flt(self._as_mpf_val(other._prec), other._mpf_))
 635             return bool(self.p * other.q < self.q * other.p)
 636         return Basic.__lt__(self, other)
 637
 638     def __le__(self, other):
 639         try:
 640             other = _sympify(other)
 641         except SympifyError:
 642             return False    # sympy > other  -->  not <=
 643         if isinstance(other, NumberSymbol):
 644             return other.__gt__(self)
 645         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 646         if isinstance(other, Number):
 647             if isinstance(other, Real):
 648                 return bool(mlib.fle(self._as_mpf_val(other._prec), other._mpf_))
 649             return bool(self.p * other.q <= self.q * other.p)
 650         return Basic.__le__(self, other)
 651
 652     def factors(self):
 653         f = factor_trial_division(self.p).copy()
 654         for p,e in factor_trial_division(self.q).items():
 655             try: f[p] += -e
 656             except KeyError: f[p] = -e
 657
 658         if len(f)>1 and f.has_key(1): del f[1]
 659         return f
 660
 661     def as_numer_denom(self):
 662         return Integer(self.p), Integer(self.q)
 663
 664     def _sage_(self):
 665         import sage.all as sage
 666         #XXX: fixme, this should work:
 667         #return sage.Integer(self[0])/sage.Integer(self[1])
 668         return sage.Integer(self.p)/sage.Integer(self.q)
 669
 670 # int -> Integer
 671 _intcache = {}
 672
 673
 674 # TODO move this tracing facility to  sympy/core/trace.py  ?
 675 def _intcache_printinfo():
 676     ints = sorted(_intcache.keys())
 677     nhit = _intcache_hits
 678     nmiss= _intcache_misses
 679
 680     if nhit == 0 and nmiss == 0:
 681         print
 682         print 'Integer cache statistic was not collected'
 683         return
 684
 685     miss_ratio = float(nmiss) / (nhit+nmiss)
 686
 687     print
 688     print 'Integer cache statistic'
 689     print '-----------------------'
 690     print
 691     print '#items: %i' % len(ints)
 692     print
 693     print ' #hit   #miss               #total'
 694     print
 695     print '%5i   %5i (%7.5f %%)   %5i'    % (nhit, nmiss, miss_ratio*100, nhit+nmiss)
 696     print
 697     print ints
 698
 699 _intcache_hits   = 0
 700 _intcache_misses = 0
 701
 702 def int_trace(f):
 703     import os
 704     if os.getenv('SYMPY_TRACE_INT', 'no').lower() != 'yes':
 705         return f
 706
 707     def Integer_tracer(cls, i):
 708         global _intcache_hits, _intcache_misses
 709
 710         try:
 711             _intcache_hits += 1
 712             return _intcache[i]
 713         except KeyError:
 714             _intcache_hits   -= 1
 715             _intcache_misses += 1
 716
 717             return f(cls, i)
 718
 719
 720     # also we want to hook our _intcache_printinfo into sys.atexit
 721     import atexit
 722     atexit.register(_intcache_printinfo)
 723
 724     return Integer_tracer
 725
 726
 727
 728
 729 class Integer(Rational):
 730
 731     q = 1
 732     is_integer = True
 733
 734     is_Integer = True
 735
 736     __slots__ = ['p']
 737
 738     def _as_mpf_val(self, prec):
 739         return mlib.from_int(self.p)
 740
 741     # TODO caching with decorator, but not to degrade performance
 742     @int_trace
 743     def __new__(cls, i):
 744         try:
 745             return _intcache[i]
 746         except KeyError:
 747             # The most often situation is when Integers are created from Python
 748             # int or long
 749             if isinstance(i, (int, long)):
 750                 obj = Basic.__new__(cls)
 751                 obj.p = i
 752                 _intcache[i] = obj
 753                 return obj
 754
 755             # Also, we seldomly need the following to work:
 756             # UC: Integer(Integer(4))   <-- sympify('4')
 757             elif i.is_Integer:
 758                 return i
 759
 760             else:
 761                 raise ValueError('invalid argument for Integer: %r' % (i,))
 762
 763     # Arithmetic operations are here for efficiency
 764     def __int__(self):
 765         return self.p
 766
 767     def __neg__(self):
 768         return Integer(-self.p)
 769
 770     def __abs__(self):
 771         if self.p >= 0:
 772             return self
 773         else:
 774             return Integer(-self.p)
 775
 776     # TODO make it decorator + bytecodehacks?
 777     def __add__(a, b):
 778         if type(b) is int:
 779             return Integer(a.p + b)
 780         elif isinstance(b, Integer):
 781             return Integer(a.p + b.p)
 782         return Rational.__add__(a, b)   # a,b -not- b,a
 783
 784     def __radd__(a, b):
 785         if type(b) is int:
 786             return Integer(b + a.p)
 787         elif isinstance(b, Integer):
 788             return Integer(b.p + a.p)
 789         return Rational.__add__(a, b)
 790
 791     def __sub__(a, b):
 792         if type(b) is int:
 793             return Integer(a.p - b)
 794         elif isinstance(b, Integer):
 795             return Integer(a.p - b.p)
 796         return Rational.__sub__(a, b)
 797
 798     def __rsub__(a, b):
 799         if type(b) is int:
 800             return Integer(b - a.p)
 801         elif isinstance(b, Integer):
 802             return Integer(b.p - a.p)
 803         return Rational.__rsub__(a, b)
 804
 805     def __mul__(a, b):
 806         if type(b) is int:
 807             return Integer(a.p * b)
 808         elif isinstance(b, Integer):
 809             return Integer(a.p * b.p)
 810         return Rational.__mul__(a, b)
 811
 812     def __rmul__(a, b):
 813         if type(b) is int:
 814             return Integer(b * a.p)
 815         elif isinstance(b, Integer):
 816             return Integer(b.p * a.p)
 817         return Rational.__mul__(a, b)
 818
 819     # XXX __pow__ ?
 820
 821     # XXX do we need to define __cmp__ ?
 822 #   def __cmp__(a, b):
 823
 824     def __eq__(a, b):
 825         if type(b) is int:
 826             return (a.p == b)
 827         elif isinstance(b, Integer):
 828             return (a.p == b.p)
 829         return Rational.__eq__(a, b)
 830
 831     def __ne__(a, b):
 832         if type(b) is int:
 833             return (a.p != b)
 834         elif isinstance(b, Integer):
 835             return (a.p != b.p)
 836         return Rational.__ne__(a, b)
 837
 838     def __gt__(a, b):
 839         if type(b) is int:
 840             return (a.p >  b)
 841         elif isinstance(b, Integer):
 842             return (a.p >  b.p)
 843         return Rational.__gt__(a, b)
 844
 845     def __lt__(a, b):
 846         if type(b) is int:
 847             return (a.p <  b)
 848         elif isinstance(b, Integer):
 849             return (a.p <  b.p)
 850         return Rational.__lt__(a, b)
 851
 852     def __ge__(a, b):
 853         if type(b) is int:
 854             return (a.p >= b)
 855         elif isinstance(b, Integer):
 856             return (a.p >= b.p)
 857         return Rational.__ge__(a, b)
 858
 859     def __le__(a, b):
 860         if type(b) is int:
 861             return (a.p <= b)
 862         elif isinstance(b, Integer):
 863             return (a.p <= b.p)
 864         return Rational.__le__(a, b)
 865
 866     ########################################
 867
 868     def _eval_is_odd(self):
 869         return bool(self.p % 2)
 870
 871     def _eval_power(b, e):
 872         if isinstance(e, Number):
 873             if e is S.NaN: return S.NaN
 874             if isinstance(e, Real):
 875                 return b._eval_evalf(e._prec) ** e
 876             if e.is_negative:
 877                 # (3/4)**-2 -> (4/3)**2
 878                 ne = -e
 879                 if ne is S.One:
 880                     return Rational(1, b.p)
 881                 return Rational(1, b.p) ** ne
 882             if e is S.Infinity:
 883                 if b.p > 1:
 884                     # (3)**oo -> oo
 885                     return S.Infinity
 886                 if b.p < -1:
 887                     # (-3)**oo -> oo + I*oo
 888                     return S.Infinity + S.Infinity * S.ImaginaryUnit
 889                 return S.Zero
 890             if isinstance(e, Integer):
 891                 # (4/3)**2 -> 4**2 / 3**2
 892                 return Integer(b.p ** e.p)
 893             if isinstance(e, Rational):
 894                 if b == -1:  # any one has tested this ???
 895                     # calculate the roots of -1
 896                     if e.q.is_odd:
 897                         return -1
 898                     r = cos(pi/e.q) + S.ImaginaryUnit*sin(pi/e.q)
 899                     return r**e.p
 900                 if b >= 0:
 901                     x, xexact = integer_nthroot(b.p, e.q)
 902                     if xexact:
 903                         res = Integer(x ** abs(e.p))
 904                         if e >= 0:
 905                             return res
 906                         else:
 907                             return 1/res
 908                     else:
 909                         if b > 2**32: #Prevent from factorizing too big integers:
 910                             for i in xrange(2, e.q/2 + 1): #OLD CODE
 911                                 if e.q % i == 0:
 912                                     x, xexact = integer_nthroot(b.p, i)
 913                                     if xexact:
 914                                         return Integer(x)**(e * i)
 915                             # Try to get some part of the base out, if exponent > 1
 916                             if e.p > e.q:
 917                                 i = e.p / e.q
 918                                 r = e.p % e.q
 919                                 return b**i * b**Rational(r, e.q)
 920                             return
 921
 922
 923                         dict = b.factors()
 924
 925                         out_int = 1
 926                         sqr_int = 1
 927                         sqr_gcd = 0
 928
 929                         sqr_dict = {}
 930
 931                         for prime,exponent in dict.iteritems():
 932                             exponent *= e.p
 933                             div_e = exponent / e.q
 934                             div_m = exponent % e.q
 935
 936                             if div_e > 0:
 937                                 out_int *= prime**div_e
 938                             if div_m > 0:
 939                                 sqr_dict[prime] = div_m
 940
 941                         for p,ex in sqr_dict.iteritems():
 942                             if sqr_gcd == 0:
 943                                 sqr_gcd = ex
 944                             else:
 945                                 sqr_gcd = igcd(sqr_gcd, ex)
 946
 947                         for k,v in sqr_dict.iteritems():
 948                             sqr_int *= k**(v/sqr_gcd)
 949
 950                         if sqr_int == b.p and out_int == 1:
 951                             return None
 952
 953                         return out_int * Pow(sqr_int , Rational(sqr_gcd, e.q))
 954                 else:
 955                     if e.q == 2:
 956                         return S.ImaginaryUnit ** e.p * (-b)**e
 957                     else:
 958                         return None
 959
 960         c,t = b.as_coeff_terms()
 961         if e.is_even and isinstance(c, Number) and c < 0:
 962             return (-c * Mul(*t)) ** e
 963
 964         return
 965
 966     def _eval_is_prime(self):
 967         if self.p < 0:
 968             return False
 969
 970     def as_numer_denom(self):
 971         return self, S.One
 972
 973     def __floordiv__(self, other):
 974         return Integer(self.p // Integer(other).p)
 975
 976     def __rfloordiv__(self, other):
 977         return Integer(Integer(other).p // self.p)
 978
 979 class Zero(Integer):
 980     __metaclass__ = SingletonMeta
 981
 982     p = 0
 983     q = 1
 984     is_positive = False
 985     is_negative = False
 986     is_finite = False
 987     is_zero = True
 988     is_prime = False
 989     is_composite = False
 990
 991     __slots__ = []
 992
 993     @staticmethod
 994     def __abs__():
 995         return S.Zero
 996
 997     @staticmethod
 998     def __neg__():
 999         return S.Zero
1000
1001     def _eval_power(b, e):
1002         if e.is_negative:
1003             return S.Infinity
1004         if e.is_positive:
1005             return b
1006         d = e.evalf()
1007         if isinstance(d, Number):
1008             if d.is_negative:
1009                 return S.Infinity
1010             return b
1011         coeff, terms = e.as_coeff_terms()
1012         if coeff.is_negative:
1013             return S.Infinity ** Mul(*terms)
1014         if coeff is not S.One:
1015             return b ** Mul(*terms)
1016
1017     def _eval_order(self, *symbols):
1018         # Order(0,x) -> 0
1019         return self
1020
1021 class One(Integer):
1022     __metaclass__ = SingletonMeta
1023
1024     p = 1
1025     q = 1
1026
1027     is_prime = True
1028
1029     __slots__ = []
1030
1031     def _eval_evalf(self, prec):
1032         return self
1033
1034     @staticmethod
1035     def __abs__():
1036         return S.One
1037
1038     @staticmethod
1039     def __neg__():
1040         return S.NegativeOne
1041
1042     def _eval_power(b, e):
1043         return b
1044
1045     def _eval_order(self, *symbols):
1046         return
1047
1048     @staticmethod
1049     def factors():
1050         return {1: 1}
1051
1052 class NegativeOne(Integer):
1053     __metaclass__ = SingletonMeta
1054
1055     p = -1
1056     q = 1
1057
1058     __slots__ = []
1059
1060     def _eval_evalf(self, prec):
1061         return self
1062
1063     @staticmethod
1064     def __abs__():
1065         return S.One
1066
1067     @staticmethod
1068     def __neg__():
1069         return S.One
1070
1071     def _eval_power(b, e):
1072         if e.is_odd: return S.NegativeOne
1073         if e.is_even: return S.One
1074         if isinstance(e, Number):
1075             if isinstance(e, Real):
1076                 return Real(-1.0) ** e
1077             if e is S.NaN:
1078                 return S.NaN
1079             if e is S.Infinity  or  e is S.NegativeInfinity:
1080                 return S.NaN
1081             if e is S.Half:
1082                 return S.ImaginaryUnit
1083             if isinstance(e, Rational):
1084                 if e.q == 2:
1085                     return S.ImaginaryUnit ** Integer(e.p)
1086                 q = int(e)
1087                 if q:
1088                     q = Integer(q)
1089                     return b ** q * b ** (e - q)
1090         return
1091
1092 class Half(Rational):
1093     __metaclass__ = SingletonMeta
1094
1095     p = 1
1096     q = 2
1097
1098     __slots__ = []
1099
1100     @staticmethod
1101     def __abs__():
1102         return S.Half
1103
1104
1105 class Infinity(Rational):
1106     __metaclass__ = SingletonMeta
1107
1108     p = 1
1109     q = 0
1110
1111     __slots__ = []
1112
1113     is_commutative = True
1114     is_positive = True
1115     is_bounded = False
1116     is_finite   = False
1117     is_infinitesimal = False
1118     is_integer  = None
1119     is_rational = None
1120     is_odd = None
1121
1122     @staticmethod
1123     def __abs__():
1124         return S.Infinity
1125
1126     @staticmethod
1127     def __neg__():
1128         return S.NegativeInfinity
1129
1130     def _eval_power(b, e):
1131         """
1132         e is symbolic object but not equal to 0, 1
1133
1134         oo ** nan -> nan
1135         oo ** (-p) -> 0, p is number, oo
1136         """
1137         if e.is_positive:
1138             return S.Infinity
1139         if e.is_negative:
1140             return S.Zero
1141         if isinstance(e, Number):
1142             if e is S.NaN:
1143                 return S.NaN
1144         d = e.evalf()
1145         if isinstance(d, Number):
1146             return b ** d
1147         return
1148
1149     def _as_mpf_val(self, prec):
1150         return mlib.finf
1151
1152     def _sage_(self):
1153         import sage.all as sage
1154         return sage.oo
1155
1156 class NegativeInfinity(Rational):
1157     __metaclass__ = SingletonMeta
1158
1159     p = -1
1160     q = 0
1161
1162     __slots__ = []
1163
1164     is_commutative = True
1165     is_real = True
1166     is_positive = False
1167     is_bounded = False
1168     is_finite = False
1169     is_infinitesimal = False
1170     is_integer  = None
1171     is_rational = None
1172
1173     @staticmethod
1174     def __abs__():
1175         return S.Infinity
1176
1177     @staticmethod
1178     def __neg__():
1179         return S.Infinity
1180
1181     def _eval_power(b, e):
1182         """
1183         e is symbolic object but not equal to 0, 1
1184
1185         (-oo) ** nan -> nan
1186         (-oo) ** oo  -> nan
1187         (-oo) ** (-oo) -> nan
1188         (-oo) ** e -> oo, e is positive even integer
1189         (-oo) ** o -> -oo, o is positive odd integer
1190
1191         """
1192         if isinstance(e, Number):
1193             if (e is S.NaN)  or  (e is S.Infinity)  or  (e is S.NegativeInfinity):
1194                 return S.NaN
1195             if isinstance(e, Integer):
1196                 if e.is_positive:
1197                     if e.is_odd:
1198                         return S.NegativeInfinity
1199                     return S.Infinity
1200             return S.NegativeOne**e * S.Infinity ** e
1201         return
1202
1203     def _as_mpf_val(self, prec):
1204         return mlib.fninf
1205
1206 class NaN(Rational):
1207     __metaclass__ = SingletonMeta
1208
1209     p = 0
1210     q = 0
1211
1212     is_commutative = True
1213     is_real = None
1214     is_rational = None
1215     is_integer  = None
1216     is_comparable = False
1217     is_finite   = None
1218     is_bounded = None
1219     #is_unbounded = False
1220     is_zero     = None
1221     is_prime    = None
1222     is_positive = None
1223
1224     __slots__ = []
1225
1226     def _as_mpf_val(self, prec):
1227         return mlib.fnan
1228
1229     def _eval_power(b, e):
1230         if e is S.Zero:
1231             return S.One
1232         return b
1233
1234     def _sage_(self):
1235         import sage.all as sage
1236         return sage.NaN
1237
1238 class ComplexInfinity(Atom):
1239     __metaclass__ = SingletonMeta
1240
1241     is_commutative = True
1242     is_comparable = None
1243     is_bounded = False
1244     is_real = None
1245
1246     __slots__ = []
1247
1248     @staticmethod
1249     def __abs__():
1250         return S.Infinity
1251
1252     @staticmethod
1253     def __neg__():
1254         return S.ComplexInfinity
1255
1256     def _eval_power(b, e):
1257         if e is S.ComplexInfinity:
1258             return S.NaN
1259
1260         if isinstance(e, Number):
1261             if e is S.Zero:
1262                 return S.NaN
1263             else:
1264                 if e.is_positive:
1265                     return S.ComplexInfinity
1266                 else:
1267                     return S.Zero
1268
1269 class NumberSymbol(Atom):
1270     __metaclass__ = SingletonMeta
1271
1272     is_commutative = True
1273     is_comparable = True
1274     is_bounded = True
1275     is_finite = True
1276
1277     __slots__ = []
1278
1279     is_NumberSymbol = True
1280
1281     def approximation(self, number_cls):
1282         """ Return an interval with number_cls endpoints
1283         that contains the value of NumberSymbol.
1284         If not implemented, then return None.
1285         """
1286
1287     def _eval_evalf(self, prec):
1288         return Real._new(self._as_mpf_val(prec), prec)
1289
1290     def _eval_derivative(self, s):
1291         return S.Zero
1292
1293     def __eq__(self, other):
1294         try:
1295             other = _sympify(other)
1296         except SympifyError:
1297             return False    # sympy != other  -->  not ==
1298         if self is other: return True
1299         if isinstance(other, Number) and self.is_irrational: return False
1300
1301         return False    # NumberSymbol != non-(Number|self)
1302
1303     def __ne__(self, other):
1304         try:
1305             other = _sympify(other)
1306         except SympifyError:
1307             return True     # sympy != other
1308         if self is other: return False
1309         if isinstance(other, Number) and self.is_irrational: return True
1310
1311         return True     # NumberSymbol != non(Number|self)
1312
1313     def __lt__(self, other):
1314         try:
1315             other = _sympify(other)
1316         except SympifyError:
1317             return False    # sympy > other  --> not <
1318         if self is other: return False
1319         if isinstance(other, Number):
1320             approx = self.approximation_interval(other.__class__)
1321             if approx is not None:
1322                 l,u = approx
1323                 if other < l: return False
1324                 if other > u: return True
1325             return self.evalf()<other
1326         if other.is_comparable:
1327             other = other.evalf()
1328             return self.evalf()<other
1329         return Basic.__lt__(self, other)
1330     def __le__(self, other):
1331         try:
1332             other = _sympify(other)
1333         except SympifyError:
1334             return False    # sympy > other  --> not <=
1335         if self is other: return True
1336         if other.is_comparable: other = other.evalf()
1337         if isinstance(other, Number):
1338             return self.evalf()<=other
1339         return Basic.__le__(self, other)
1340     def __gt__(self, other):
1341         return (-self) < (-other)
1342     def __ge__(self, other):
1343         return (-self) <= (-other)
1344
1345
1346 class Exp1(NumberSymbol):
1347
1348     is_real = True
1349     is_positive = True
1350     is_negative = False # XXX Forces is_negative/is_nonnegative
1351     is_irrational = True
1352
1353     __slots__ = []
1354
1355     @staticmethod
1356     def __abs__():
1357         return S.Exp1
1358
1359     def _as_mpf_val(self, prec):
1360         return mlib.fe(prec)
1361
1362     def approximation_interval(self, number_cls):
1363         if issubclass(number_cls,Integer):
1364             return (Integer(2),Integer(3))
1365         elif issubclass(number_cls,Rational):
1366             pass
1367
1368     def _eval_power(self, exp):
1369         return C.exp(exp)
1370
1371     def _sage_(self):
1372         import sage.all as sage
1373         return sage.e
1374
1375 class Pi(NumberSymbol):
1376
1377     is_real = True
1378     is_positive = True
1379     is_negative = False
1380     is_irrational = True
1381
1382     __slots__ = []
1383
1384     @staticmethod
1385     def __abs__():
1386         return S.Pi
1387
1388     def _as_mpf_val(self, prec):
1389         return mlib.fpi(prec)
1390
1391     def approximation_interval(self, number_cls):
1392         if issubclass(number_cls, Integer):
1393             return (Integer(3), Integer(4))
1394         elif issubclass(number_cls, Rational):
1395             return (Rational(223,71), Rational(22,7))
1396
1397     def _sage_(self):
1398         import sage.all as sage
1399         return sage.pi
1400
1401 class GoldenRatio(NumberSymbol):
1402
1403     is_real = True
1404     is_positive = True
1405     is_negative = False
1406     is_irrational = True
1407
1408     __slots__ = []
1409
1410     def _as_mpf_val(self, prec):
1411         return mlib.from_man_exp(mpmath.specfun.phi_fixed(prec+10), -prec-10)
1412
1413     def _eval_expand_func(self, *args):
1414         return S.Half + S.Half*S.Sqrt(5)
1415
1416     def approximation_interval(self, number_cls):
1417         if issubclass(number_cls, Integer):
1418             return (S.One, Rational(2))
1419         elif issubclass(number_cls, Rational):
1420             pass
1421
1422     def _sage_(self):
1423         import sage.all as sage
1424         return sage.golden_ratio
1425
1426 class EulerGamma(NumberSymbol):
1427
1428     is_real = True
1429     is_positive = True
1430     is_negative = False
1431     is_irrational = None
1432
1433     __slots__ = []
1434
1435     def _as_mpf_val(self, prec):
1436         return mlib.from_man_exp(mpmath.specfun.euler_fixed(prec+10), -prec-10)
1437
1438     def approximation_interval(self, number_cls):
1439         if issubclass(number_cls, Integer):
1440             return (S.Zero, S.One)
1441         elif issubclass(number_cls, Rational):
1442             return (S.Half, Rational(3, 5))
1443
1444     def _sage_(self):
1445         import sage.all as sage
1446         return sage.euler_gamma
1447
1448 class Catalan(NumberSymbol):
1449
1450     is_real = True
1451     is_positive = True
1452     is_negative = False
1453     is_irrational = None
1454
1455     __slots__ = []
1456
1457     def _as_mpf_val(self, prec):
1458         return mlib.from_man_exp(mpmath.specfun.catalan_fixed(prec+10), -prec-10)
1459
1460     def approximation_interval(self, number_cls):
1461         if issubclass(number_cls, Integer):
1462             return (S.Zero, S.One)
1463         elif issubclass(number_cls, Rational):
1464             return (Rational(9, 10), S.One)
1465
1466     def _sage_(self):
1467         import sage.all as sage
1468         return sage.catalan
1469
1470 class ImaginaryUnit(Atom):
1471     __metaclass__ = SingletonMeta
1472
1473     is_commutative = True
1474     is_imaginary = True
1475     is_bounded = True
1476     is_finite = True
1477
1478     __slots__ = []
1479
1480     @staticmethod
1481     def __abs__():
1482         return S.One
1483
1484     def _eval_evalf(self, prec):
1485         return self
1486
1487     def _eval_conjugate(self):
1488         return -S.ImaginaryUnit
1489
1490     def _eval_derivative(self, s):
1491         return S.Zero
1492
1493     def _eval_power(b, e):
1494         """
1495         b is I = sqrt(-1)
1496         e is symbolic object but not equal to 0, 1
1497
1498         I ** r -> (-1)**(r/2) -> exp(r/2 * Pi * I) -> sin(Pi*r/2) + cos(Pi*r/2) * I, r is decimal
1499         I ** 0 mod 4 -> 1
1500         I ** 1 mod 4 -> I
1501         I ** 2 mod 4 -> -1
1502         I ** 3 mod 4 -> -I
1503         """
1504
1505
1506         if isinstance(e, Number):
1507             #if isinstance(e, Decimal):
1508             #    a = decimal_math.pi() * exponent.num / 2
1509             #    return Decimal(decimal_math.sin(a) + decimal_math.cos(a) * ImaginaryUnit())
1510             if isinstance(e, Integer):
1511                 e = e.p % 4
1512                 if e==0: return S.One
1513                 if e==1: return S.ImaginaryUnit
1514                 if e==2: return -S.One
1515                 return -S.ImaginaryUnit
1516             return (-S.One) ** (e * S.Half)
1517         return
1518
1519     def as_base_exp(self):
1520         return -S.One, S.Half
1521
1522     def _sage_(self):
1523         import sage.all as sage
1524         return sage.I
1525
1526 # /cyclic/
1527 import basic as _
1528 _.Number    = Number
1529 _.Integer   = Integer
1530 _.Rational  = Rational
1531 _.Real      = Real
1532 del _
1533
1534 import add as _
1535 _.Number    = Number
1536 del _
1537
1538 import mul as _
1539 _.Number    = Number
1540 _.Integer   = Integer
1541 _.Real      = Real
1542 del _
1543
1544 import power as _
1545 _.Number    = Number
1546 _.Rational  = Rational
1547 _.Integer   = Integer
1548 del _
1549
1550 import sympify as _
1551 _.Integer   = Integer
1552 _.Real      = Real
1553 del _
1554
1555 # ----
1556
1557 _intcache[0] = S.Zero
1558 _intcache[1] = S.One
1559 _intcache[-1]= S.NegativeOne
1560
1561 Basic.singleton['E'] = Exp1
1562 Basic.singleton['pi'] = Pi
1563 Basic.singleton['I'] = ImaginaryUnit
1564 Basic.singleton['oo'] = Infinity
1565 Basic.singleton['nan'] = NaN
1566
1567 Basic.singleton['zoo'] = ComplexInfinity
1568
1569 Basic.singleton['GoldenRatio'] = GoldenRatio
1570 Basic.singleton['EulerGamma'] = EulerGamma
1571 Basic.singleton['Catalan'] = Catalan