sympy/core/numbers.py

   1 import math
   2 import decimal
   3 import decimal_math
   4 from basic import Basic, Atom, Singleton, S, Memoizer, MemoizerArg
   5 from methods import NoRelMeths, RelMeths, ArithMeths
   6 from power import integer_nthroot
   7
   8 @Memoizer((int, long), (int, long))
   9 def gcd(a, b):
  10     """
  11     Returns the Greatest Common Divisor, implementing Euclid's algorithm.
  12     """
  13     while a:
  14         a, b = b%a, a
  15     return b
  16
  17 @Memoizer((int, long), return_value_converter = lambda d: d.copy())
  18 def factor_trial_division(n):
  19     """
  20     Factor any integer into a product of primes, 0, 1, and -1.
  21     Returns a dictionary {<prime: exponent>}.
  22     """
  23     if not n:
  24         return {0:1}
  25     factors = {}
  26     if n < 0:
  27         factors[-1] = 1
  28         n = -n
  29     if n==1:
  30         factors[1] = 1
  31         return factors
  32     d = 2
  33     while n % d == 0:
  34         try:
  35             factors[d] += 1
  36         except KeyError:
  37             factors[d] = 1
  38         n //= d
  39     d = 3
  40     while n > 1 and d*d <= n:
  41         if n % d:
  42             d += 2
  43         else:
  44             try:
  45                 factors[d] += 1
  46             except KeyError:
  47                 factors[d] = 1
  48             n //= d
  49     if n>1:
  50         try:
  51             factors[n] += 1
  52         except KeyError:
  53             factors[n] = 1
  54     return factors
  55
  56
  57 class Number(Atom, RelMeths, ArithMeths):
  58     """
  59     Represents any kind of number in sympy.
  60
  61     Floating point numbers are represented by the Real class.
  62     Integer numbers (of any size), together with rational numbers (again, there
  63     is no limit on their size) are represented by the Rational class.
  64
  65     If you want to represent for example 1+sqrt(2), then you need to do:
  66
  67     Rational(1) + sqrt(Rational(2))
  68     """
  69     is_commutative = True
  70     is_comparable = True
  71     is_bounded = True
  72     is_finite = True
  73
  74     def __new__(cls, *obj):
  75         if len(obj)==1: obj=obj[0]
  76         if isinstance(obj, (int, long)):
  77             return Integer(obj)
  78         if isinstance(obj,tuple) and len(obj)==2:
  79             return Rational(*obj)
  80         if isinstance(obj, (str,float,decimal.Decimal)):
  81             return Real(obj)
  82         if isinstance(obj, Number):
  83             return obj
  84         raise TypeError("expected str|int|long|float|Decimal|Number object but got %r" % (obj))
  85
  86     def _eval_evalf(self):
  87         r = self._as_decimal()
  88         return Real(r)
  89
  90     def __float__(self):
  91         return float(self._as_decimal())
  92
  93     def _as_decimal(self):
  94         raise NotImplementedError('%s needs ._as_decimal() method' % (self.__class__.__name__))
  95
  96     def _eval_derivative(self, s):
  97         return S.Zero
  98
  99     def _eval_conjugate(self):
 100         return self
 101
 102     def _eval_apply(self, a):
 103         return self*a
 104
 105     def _eval_order(self, *symbols):
 106         # Order(5, x, y) -> Order(1,x,y)
 107         return Basic.Order(S.One, *symbols)
 108
 109     def sqrt(self): return Real(decimal_math.sqrt(self._as_decimal()))
 110     def exp(self): return Real(decimal_math.exp(self._as_decimal()))
 111     def log(self): return Real(decimal_math.log(self._as_decimal()))
 112     def sin(self): return Real(decimal_math.sin(self._as_decimal()))
 113     def cos(self): return Real(decimal_math.cos(self._as_decimal()))
 114     def tan(self): return Real(decimal_math.tan(self._as_decimal()))
 115     def cot(self): return Real(decimal_math.cot(self._as_decimal()))
 116     def asin(self): return Real(decimal_math.asin(self._as_decimal()))
 117     def acos(self): return Real(decimal_math.acos(self._as_decimal()))
 118     def atan(self): return Real(decimal_math.atan(self._as_decimal()))
 119     def acot(self): return Real(decimal_math.acot(self._as_decimal()))
 120     def sinh(self): return Real(decimal_math.sinh(self._as_decimal()))
 121     def cosh(self): return Real(decimal_math.cosh(self._as_decimal()))
 122     def tanh(self): return Real(decimal_math.tanh(self._as_decimal()))
 123     def coth(self): return Real(decimal_math.coth(self._as_decimal()))
 124     def asinh(self): return Real(decimal_math.asinh(self._as_decimal()))
 125     def acosh(self): return Real(decimal_math.acosh(self._as_decimal()))
 126     def atanh(self): return Real(decimal_math.atanh(self._as_decimal()))
 127     def acoth(self): return Real(decimal_math.acoth(self._as_decimal()))
 128     def floor(self): return Real(decimal_math.floor(self._as_decimal()))
 129     def ceiling(self): return Real(decimal_math.ceiling(self._as_decimal()))
 130
 131     def __eq__(self, other):
 132         raise NotImplementedError,'%s needs .__eq__() method' % (self.__class__.__name__)
 133     def __ne__(self, other):
 134         raise NotImplementedError,'%s needs .__ne__() method' % (self.__class__.__name__)
 135     def __lt__(self, other):
 136         raise NotImplementedError,'%s needs .__lt__() method' % (self.__class__.__name__)
 137     def __le__(self, other):
 138         raise NotImplementedError,'%s needs .__le__() method' % (self.__class__.__name__)
 139
 140     def __gt__(self, other):
 141         return Basic.sympify(other).__lt__(self)
 142     def __ge__(self, other):
 143         return Basic.sympify(other).__le__(self)
 144
 145     def as_coeff_terms(self, x=None):
 146         # a -> c * t
 147         return self, []
 148
 149 decimal_to_Number_cls = {
 150     decimal.Decimal('0').as_tuple():'Zero',
 151     decimal.Decimal('-0').as_tuple():'Zero',
 152     decimal.Decimal('1').as_tuple():'One',
 153     decimal.Decimal('-1').as_tuple():'NegativeOne',
 154     decimal.Decimal('Infinity').as_tuple():'Infinity',
 155     decimal.Decimal('-Infinity').as_tuple():'NegativeInfinity',
 156     decimal.Decimal('NaN').as_tuple():'NaN',
 157     }
 158
 159 def convert_to_Decimal(num):
 160     if isinstance(num, str):
 161         num = decimal.Decimal(num)
 162     elif isinstance(num, float):
 163         num = Real.float_to_decimal(num)
 164     return +num
 165
 166 class Real(Number):
 167     """
 168     Represents a floating point number. It is capable of representing
 169     arbitrary-precision floating-point numbers
 170
 171     Usage:
 172     ======
 173         Real(3.5)   .... 3.5 (the 3.5 was converted from a python float)
 174         Real("3.0000000000000005")
 175
 176     Notes:
 177     ======
 178         - Real(x) with x being a Python int/long will return Integer(x)
 179     """
 180     is_real = True
 181     is_irrational = False
 182     is_integer = False
 183
 184     @Memoizer(type, MemoizerArg((str, int, long, float, decimal.Decimal), convert_to_Decimal))
 185     def __new__(cls, num):
 186         if isinstance(num, (int, long)):
 187             return Integer(num)
 188
 189         singleton_cls_name = decimal_to_Number_cls.get(num.as_tuple(), None)
 190         if singleton_cls_name is not None:
 191             return getattr(Basic, singleton_cls_name)()
 192         obj = Basic.__new__(cls)
 193         obj.num = num
 194         return obj
 195
 196     @staticmethod
 197     def float_to_decimal(f):
 198         "Convert a floating point number to a Decimal with no loss of information"
 199         # Transform (exactly) a float to a mantissa (0.5 <= abs(m) < 1.0) and an
 200         # exponent.  Double the mantissa until it is an integer.  Use the integer
 201         # mantissa and exponent to compute an equivalent Decimal.  If this cannot
 202         # be done exactly, then retry with more precision.
 203
 204         try:
 205             mantissa, exponent = math.frexp(f)
 206         except OverflowError:
 207             return decimal.Inf
 208
 209         while mantissa != int(mantissa):
 210             mantissa *= 2.0
 211             exponent -= 1
 212         mantissa = int(mantissa)
 213
 214         oldcontext = decimal.getcontext()
 215         decimal.setcontext(decimal.Context(traps=[decimal.Inexact]))
 216         try:
 217             while True:
 218                 try:
 219                     return mantissa * decimal.Decimal(2) ** exponent
 220                 except decimal.Inexact:
 221                     decimal.getcontext().prec += 1
 222         finally:
 223             decimal.setcontext(oldcontext)
 224
 225     def _hashable_content(self):
 226         return (self.num,)
 227
 228     def tostr(self, level=0):
 229         r = str(self.num.normalize())
 230         if self.precedence<=level:
 231             return '(%s)' % (r)
 232         return r
 233
 234     def torepr(self):
 235         return '%s(%r)' % (self.__class__.__name__, str(self.num))
 236
 237     def _eval_is_positive(self):
 238         return self.num.as_tuple()[0] == 0
 239
 240     def _eval_is_negative(self):
 241         return self.num.as_tuple()[0] == 1
 242
 243     def _eval_evalf(self):
 244         return self
 245
 246     def _as_decimal(self):
 247         return self.num
 248
 249     def __neg__(self):
 250         return Real(-self.num)
 251
 252     def __mul__(self, other):
 253         other = Basic.sympify(other)
 254         if isinstance(other, Number):
 255             return Real(self.num * other._as_decimal())
 256         return Number.__mul__(self, other)
 257
 258     def __add__(self, other):
 259         other = Basic.sympify(other)
 260         if isinstance(other, NaN) or isinstance(self, NaN):
 261             return S.NaN
 262         if isinstance(other, Number):
 263             return Real(self.num + other._as_decimal())
 264         return Number.__add__(self, other)
 265
 266     def _eval_power(b, e):
 267         """
 268         b is Real but not equal to rationals, integers, 0.5, oo, -oo, nan
 269         e is symbolic object but not equal to 0, 1
 270
 271         (-p) ** r -> exp(r * log(-p)) -> exp(r * (log(p) + I*Pi)) ->
 272                   -> p ** r * (sin(Pi*r) + cos(Pi*r) * I)
 273         """
 274         if isinstance(e, Number):
 275             if isinstance(e, Integer):
 276                 e = e.p
 277                 return Real(decimal_math.pow(b.num, e))
 278
 279             e2 = e._as_decimal()
 280             if b.is_negative and not e.is_integer:
 281                 m = decimal_math.pow(-b.num, e2)
 282                 a = decimal_math.pi() * e2
 283                 s = m * decimal_math.sin(a)
 284                 c = m * decimal_math.cos(a)
 285                 return Real(s) + Real(c) * S.ImaginaryUnit
 286             return Real(decimal_math.pow(b.num, e2))
 287         return
 288
 289     def __abs__(self):
 290         return Real(abs(self.num))
 291
 292     def __int__(self):
 293         return int(self.num)
 294
 295     def __float__(self):
 296         return float(self.num)
 297
 298     def __eq__(self, other):
 299         other = Basic.sympify(other)
 300         if isinstance(other, NumberSymbol):
 301             if other.is_irrational: return False
 302             return other.__eq__(self)
 303         if other.is_comparable: other = other.evalf()
 304         if isinstance(other, Number):
 305             return bool(self._as_decimal()==other._as_decimal())
 306         return RelMeths.__eq__(self, other)
 307     def __ne__(self, other):
 308         other = Basic.sympify(other)
 309         if isinstance(other, NumberSymbol):
 310             if other.is_irrational: return True
 311             return other.__ne__(self)
 312         if other.is_comparable: other = other.evalf()
 313         if isinstance(other, Number):
 314             return bool(self._as_decimal()!=other._as_decimal())
 315         return RelMeths.__ne__(self, other)
 316     def __lt__(self, other):
 317         other = Basic.sympify(other)
 318         if isinstance(other, NumberSymbol):
 319             return other.__ge__(self)
 320         if other.is_comparable: other = other.evalf()
 321         if isinstance(other, Number):
 322             return bool(self._as_decimal() < other._as_decimal())
 323         return RelMeths.__lt__(self, other)
 324     def __le__(self, other):
 325         other = Basic.sympify(other)
 326         if isinstance(other, NumberSymbol):
 327             return other.__gt__(self)
 328         if other.is_comparable: other = other.evalf()
 329         if isinstance(other, Number):
 330             return bool(self._as_decimal()<=other._as_decimal())
 331         return RelMeths.__le__(self, other)
 332
 333     def epsilon_eq(self, other, epsilon="10e-16"):
 334         return abs(self - other) < Basic.Real(epsilon)
 335
 336 def _parse_rational(s):
 337     """Parse rational number from string representation"""
 338     # Simple fraction
 339     if "/" in s:
 340         p, q = s.split("/")
 341         return int(p), int(q)
 342     # Recurring decimal
 343     elif "[" in s:
 344         sign = 1
 345         if s[0] == "-":
 346             sign = -1
 347             s = s[1:]
 348         s, periodic = s.split("[")
 349         periodic = periodic.rstrip("]")
 350         offset = len(s) - s.index(".") - 1
 351         n1 = int(periodic)
 352         n2 = int("9" * len(periodic))
 353         r = Rational(*_parse_rational(s)) + Rational(n1, n2*10**offset)
 354         return sign*r.p, r.q
 355     # Ordinary decimal string. Use the Decimal class's built-in parser
 356     else:
 357         sign, digits, expt = decimal.Decimal(s).as_tuple()
 358         p = (1, -1)[sign] * int("".join(str(x) for x in digits))
 359         if expt >= 0:
 360             return p*(10**expt), 1
 361         else:
 362             return p, 10**-expt
 363
 364
 365 class Rational(Number):
 366     """Represents integers and rational numbers (p/q) of any size.
 367
 368     Examples
 369     ========
 370         >>> Rational(3)
 371         3
 372         >>> Rational(1,2)
 373         1/2
 374
 375     You can create a rational from a string:
 376         >>> Rational("3/5")
 377         3/5
 378         >>> Rational("1.23")
 379         123/100
 380
 381     Use square brackets to indicate a recurring decimal:
 382         >>> Rational("0.[333]")
 383         1/3
 384         >>> Rational("1.2[05]")
 385         1193/990
 386         >>> float(Rational(1193,990))
 387         1.2050505050505051
 388     """
 389     is_real = True
 390     is_integer = False
 391     is_rational = True
 392
 393     @Memoizer(type, (int, long, str), MemoizerArg((int, long, type(None)), name="q"))
 394     def __new__(cls, p, q = None):
 395         if q is None:
 396             if isinstance(p, str):
 397                 p, q = _parse_rational(p)
 398             else:
 399                 return Integer(p)
 400         if q==0:
 401             if p==0: return S.NaN
 402             if p<0: return S.NegativeInfinity
 403             return S.Infinity
 404         if q<0:
 405             q = -q
 406             p = -p
 407         n = gcd(abs(p), q)
 408         if n>1:
 409             p /= n
 410             q /= n
 411         if q==1: return Integer(p)
 412         if p==1 and q==2: return S.Half
 413         obj = Basic.__new__(cls)
 414         obj.p = p
 415         obj.q = q
 416         return obj
 417
 418     def _hashable_content(self):
 419         return (self.p, self.q)
 420
 421     def tostr(self, level=0):
 422         if self.precedence<=level:
 423             return '(%s/%s)' % (self.p, self.q)
 424         return '%s/%s' % (self.p, self.q)
 425
 426     def torepr(self):
 427         return '%s(%r, %r)' % (self.__class__.__name__, self.p, self.q)
 428
 429     @property
 430     def precedence(self):
 431         if self.p < 0:
 432             return Basic.Add_precedence
 433         return Basic.Mul_precedence
 434
 435     def _eval_is_positive(self):
 436         return self.p > 0
 437
 438     def _eval_is_zero(self):
 439         return self.p == 0
 440
 441     def __neg__(self): return Rational(-self.p, self.q)
 442
 443     def __mul__(self, other):
 444         other = Basic.sympify(other)
 445         if isinstance(other, NaN) or isinstance(self, NaN):
 446             return S.NaN
 447         if isinstance(other, Real):
 448             return Real(self._as_decimal() * other.num)
 449         if isinstance(other, Rational):
 450             return Rational(self.p * other.p, self.q * other.q)
 451         return Number.__mul__(self, other)
 452
 453     def __add__(self, other):
 454         other = Basic.sympify(other)
 455         if isinstance(other, NaN) or isinstance(self, NaN):
 456             return S.NaN
 457         if isinstance(other, Real):
 458             return Real(self._as_decimal() + other.num)
 459         if isinstance(other, Rational):
 460             if self.is_unbounded:
 461                 if other.is_bounded:
 462                     return self
 463                 elif self==other:
 464                     return self
 465             else:
 466                 if other.is_unbounded:
 467                     return other
 468             return Rational(self.p * other.q + self.q * other.p, self.q * other.q)
 469         return Number.__add__(self, other)
 470
 471     def _eval_power(b, e):
 472         if isinstance(e, Number):
 473             if isinstance(e, NaN): return S.NaN
 474             if isinstance(e, Real):
 475                 return Real(decimal_math.pow(b._as_decimal(), e.num))
 476             if e.is_negative:
 477                 # (3/4)**-2 -> (4/3)**2
 478                 ne = -e
 479                 if isinstance(ne, One):
 480                     return Rational(b.q, b.p)
 481                 return Rational(b.q, b.p) ** ne
 482             if isinstance(e, Infinity):
 483                 if b.p > b.q:
 484                     # (3/2)**oo -> oo
 485                     return S.Infinity
 486                 if b.p < -b.q:
 487                     # (-3/2)**oo -> oo + I*oo
 488                     return S.Infinity + S.Infinity * S.ImaginaryUnit
 489                 return S.Zero
 490             if isinstance(e, Integer):
 491                 # (4/3)**2 -> 4**2 / 3**2
 492                 return Rational(b.p ** e.p, b.q ** e.p)
 493             if isinstance(e, Rational):
 494                 if b.p != 1:
 495                     # (4/3)**(5/6) -> 4**(5/6) * 3**(-5/6)
 496                     return Integer(b.p) ** e * Integer(b.q) ** (-e)
 497                 if b >= 0:
 498                     x, xexact = integer_nthroot(b.p, e.q)
 499                     y, yexact = integer_nthroot(b.q, e.q)
 500                     if xexact and yexact:
 501                         res = Rational(x ** abs(e.p), y ** abs(e.p))
 502                         if e >= 0:
 503                             return res
 504                         else:
 505                             return 1/res
 506                     # Now check also devisors of the exponents denominator
 507                     # TODO: Check if this slows down to much.
 508                     for i in xrange(2, e.q/2 + 1):
 509                         if e.q % i == 0:
 510                             x, xexact = integer_nthroot(b.p, i)
 511                             y, yexact = integer_nthroot(b.q, i)
 512                             if xexact and yexact:
 513                                 return Rational(x, y)**Rational(e.p, e.q/i)
 514                     else:
 515                         # Try to get some part of the base out, if exp > 1
 516                         if e.p > e.q:
 517                             i = e.p / e.q
 518                             r = e.p % e.q
 519                             return b**i * b**Rational(r, e.q)
 520                 else:
 521                     return (-1)**e * (-b)**e
 522
 523         c,t = b.as_coeff_terms()
 524         if e.is_even and isinstance(c, Basic.Number) and c < 0:
 525             return (-c * Basic.Mul(*t)) ** e
 526
 527         return
 528
 529     def _as_decimal(self):
 530         return decimal.Decimal(self.p) / decimal.Decimal(self.q)
 531
 532     def __abs__(self):
 533         return Rational(abs(self.p), self.q)
 534
 535     def __int__(self):
 536         return int(self.p//self.q)
 537
 538     def __eq__(self, other):
 539         other = Basic.sympify(other)
 540         if isinstance(other, NumberSymbol):
 541             if other.is_irrational: return False
 542             return other.__eq__(self)
 543         from sympy.core.function import FunctionClass
 544         if isinstance(self, Number) and isinstance(other, FunctionClass):
 545             return False
 546         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 547         if isinstance(other, Number):
 548             if isinstance(other, Real):
 549                 return bool(self._as_decimal()==other._as_decimal())
 550             return bool(self.p==other.p and self.q==other.q)
 551         return RelMeths.__eq__(self, other)
 552     def __ne__(self, other):
 553         other = Basic.sympify(other)
 554         if isinstance(other, NumberSymbol):
 555             if other.is_irrational: return True
 556             return other.__ne__(self)
 557         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 558         if isinstance(other, Number):
 559             if isinstance(other, Real):
 560                 return bool(self._as_decimal()!=other._as_decimal())
 561             return bool(self.p!=other.p or self.q!=other.q)
 562         return RelMeths.__ne__(self, other)
 563     def __lt__(self, other):
 564         other = Basic.sympify(other)
 565         if isinstance(other, NumberSymbol):
 566             return other.__ge__(self)
 567         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 568         if isinstance(other, Number):
 569             if isinstance(other, Real):
 570                 return bool(self._as_decimal() < other._as_decimal())
 571             return bool(self.p * other.q < self.q * other.p)
 572         return RelMeths.__lt__(self, other)
 573     def __le__(self, other):
 574         other = Basic.sympify(other)
 575         if isinstance(other, NumberSymbol):
 576             return other.__gt__(self)
 577         if other.is_comparable and not isinstance(other, Rational): other = other.evalf()
 578         if isinstance(other, Number):
 579             if isinstance(other, Real):
 580                 return bool(self._as_decimal()<=other._as_decimal())
 581             return bool(self.p * other.q <= self.q * other.p)
 582         return RelMeths.__le__(self, other)
 583
 584     def factors(self):
 585         f = factor_trial_division(self.p).copy()
 586         for p,e in factor_trial_division(self.q).items():
 587             try: f[p] += -e
 588             except KeyError: f[p] = -e
 589         fi = {}
 590         for p,e in f.items():
 591             if e==0:
 592                 del f[p]
 593             else:
 594                 try: fi[e] *= p
 595                 except KeyError: fi[e] = p
 596         f = {}
 597         for e,p in fi.items():
 598             f[p] = e
 599         if len(f)>1 and f.has_key(1): del f[1]
 600         return f
 601
 602     def as_numer_denom(self):
 603         return Integer(self.p), Integer(self.q)
 604
 605 class Integer(Rational):
 606
 607     q = 1
 608     is_integer = True
 609
 610     @Memoizer(type, (int, long))
 611     def __new__(cls, i):
 612         if isinstance(i, Integer):
 613             return i
 614         if i==0: return S.Zero
 615         if i==1: return S.One
 616         if i==-1: return S.NegativeOne
 617         obj = Basic.__new__(cls)
 618         obj.p = i
 619         return obj
 620
 621     def _eval_is_odd(self):
 622         return bool(self.p % 2)
 623
 624     @property
 625     def precedence(self):
 626         if self.p < 0:
 627             return 40 # same as Add
 628         return Atom.precedence
 629
 630     def tostr(self, level=0):
 631         if self.precedence<=level:
 632             return '(%s)' % (self.p)
 633         return str(self.p)
 634
 635     def torepr(self):
 636         return '%s(%r)' % (self.__class__.__name__, self.p)
 637
 638     def _eval_evalf(self):
 639         return self
 640
 641     def _eval_power(b, e):
 642         if isinstance(e, Number):
 643             if isinstance(e, NaN): return S.NaN
 644             if isinstance(e, Real):
 645                 return Real(decimal_math.pow(b._as_decimal(), e.num))
 646             if e.is_negative:
 647                 # (3/4)**-2 -> (4/3)**2
 648                 ne = -e
 649                 if isinstance(ne, One):
 650                     return Rational(1, b.p)
 651                 return Rational(1, b.p) ** ne
 652             if isinstance(e, Infinity):
 653                 if b.p > 1:
 654                     # (3)**oo -> oo
 655                     return S.Infinity
 656                 if b.p < -1:
 657                     # (-3)**oo -> oo + I*oo
 658                     return S.Infinity + S.Infinity * S.ImaginaryUnit
 659                 return S.Zero
 660             if isinstance(e, Integer):
 661                 # (4/3)**2 -> 4**2 / 3**2
 662                 return Integer(b.p ** e.p)
 663             if isinstance(e, Rational):
 664                 if b == -1:  # any one has tested this ???
 665                     # calculate the roots of -1
 666                     if e.q.is_odd:
 667                         return -1
 668                     r = cos(pi/e.q) + S.ImaginaryUnit*sin(pi/e.q)
 669                     return r**e.p
 670                 if b >= 0:
 671                     x, xexact = integer_nthroot(b.p, e.q)
 672                     if xexact:
 673                         res = Integer(x ** abs(e.p))
 674                         if e >= 0:
 675                             return res
 676                         else:
 677                             return 1/res
 678                     # Now check also devisors of the exponents denominator
 679                     # TODO: Check if this slows down to much.
 680                     for i in xrange(2, e.q/2 + 1):
 681                         if e.q % i == 0:
 682                             x, xexact = integer_nthroot(b.p, i)
 683                             if xexact:
 684                                 return Integer(x)**(e * i)
 685                     # Try to get some part of the base out, if exponent > 1
 686                     if e.p > e.q:
 687                         i = e.p / e.q
 688                         r = e.p % e.q
 689                         return b**i * b**Rational(r, e.q)
 690
 691                 else:
 692                     return (-1)**e * (-b)**e
 693
 694         c,t = b.as_coeff_terms()
 695         if e.is_even and isinstance(c, Basic.Number) and c < 0:
 696             return (-c * Basic.Mul(*t)) ** e
 697
 698         return
 699
 700     def _eval_is_prime(self):
 701         if self.p < 0:
 702             return False
 703
 704     def as_numer_denom(self):
 705         return self, S.One
 706
 707     def __floordiv__(self, other):
 708         return Integer(self.p // Integer(other).p)
 709
 710     def __rfloordiv__(self, other):
 711         return Integer(Integer(other).p // self.p)
 712
 713 class Zero(Singleton, Integer):
 714
 715     p = 0
 716     q = 1
 717     is_positive = False
 718     is_negative = False
 719     is_finite = False
 720     is_zero = True
 721     is_prime = False
 722     is_composite = True
 723
 724     def _eval_power(b, e):
 725         if e.is_negative:
 726             return S.Infinity
 727         if e.is_positive:
 728             return b
 729         d = e.evalf()
 730         if isinstance(d, Number):
 731             if d.is_negative:
 732                 return S.Infinity
 733             return b
 734         coeff, terms = e.as_coeff_terms()
 735         if coeff.is_negative:
 736             return S.Infinity ** Basic.Mul(*terms)
 737         if not isinstance(coeff, Basic.One):
 738             return b ** Basic.Mul(*terms)
 739
 740     def _eval_order(self, *symbols):
 741         # Order(0,x) -> 0
 742         return self
 743
 744 class One(Singleton, Integer):
 745
 746     p = 1
 747     q = 1
 748
 749     is_prime = True
 750
 751     def _eval_power(b, e):
 752         return b
 753
 754     def _eval_order(self, *symbols):
 755         return
 756
 757 class NegativeOne(Singleton, Integer):
 758
 759     p = -1
 760     q = 1
 761
 762     def _eval_power(b, e):
 763         if e.is_odd: return S.NegativeOne
 764         if e.is_even: return S.One
 765         if isinstance(e, Number):
 766             if isinstance(e, Real):
 767                 a = e.num * decimal_math.pi()
 768                 s = decimal_math.sin(a)
 769                 c = decimal_math.cos(a)
 770                 return Real(s) + Real(c) * S.ImaginaryUnit
 771             if isinstance(e, NaN):
 772                 return S.NaN
 773             if isinstance(e, (Infinity, NegativeInfinity)):
 774                 return S.NaN
 775             if isinstance(e, Half):
 776                 return S.ImaginaryUnit
 777             if isinstance(e, Rational):
 778                 if e.q == 2:
 779                     return S.ImaginaryUnit ** Integer(e.p)
 780                 q = int(e)
 781                 if q:
 782                     q = Integer(q)
 783                     return b ** q * b ** (e - q)
 784         return
 785
 786 class Half(Singleton, Rational):
 787
 788     p = 1
 789     q = 2
 790
 791 class Infinity(Singleton, Rational):
 792
 793     p = 1
 794     q = 0
 795
 796     is_commutative = True
 797     is_positive = True
 798     is_bounded = False
 799     is_finite = None
 800     is_odd = None
 801
 802     def tostr(self, level=0):
 803         return 'oo'
 804
 805     def _eval_power(b, e):
 806         """
 807         e is symbolic object but not equal to 0, 1
 808
 809         oo ** nan -> nan
 810         oo ** (-p) -> 0, p is number, oo
 811         """
 812         if e.is_positive:
 813             return S.Infinity
 814         if e.is_negative:
 815             return S.Zero
 816         if isinstance(e, Number):
 817             if isinstance(e, NaN):
 818                 return S.NaN
 819         d = e.evalf()
 820         if isinstance(d, Number):
 821             return b ** d
 822         return
 823
 824     def _as_decimal(self):
 825         return decimal.Decimal('Infinity')
 826
 827 class NegativeInfinity(Singleton, Rational):
 828
 829     p = -1
 830     q = 0
 831
 832     is_commutative = True
 833     is_real = True
 834     is_positive = False
 835     is_bounded = False
 836     is_finite = False
 837
 838     precedence = 40 # same as Add
 839
 840     def tostr(self, level=0):
 841         return '-oo'
 842
 843     def _eval_power(b, e):
 844         """
 845         e is symbolic object but not equal to 0, 1
 846
 847         (-oo) ** nan -> nan
 848         (-oo) ** oo  -> nan
 849         (-oo) ** (-oo) -> nan
 850         (-oo) ** e -> oo, e is positive even integer
 851         (-oo) ** o -> -oo, o is positive odd integer
 852
 853         """
 854         if isinstance(e, Number):
 855             if isinstance(e, (NaN, Infinity, NegativeInfinity)):
 856                 return S.NaN
 857             if isinstance(e, Integer):
 858                 if e.is_positive:
 859                     if e.is_odd:
 860                         return S.NegativeInfinity
 861                     return S.Infinity
 862             return S.NegativeOne**e * S.Infinity ** e
 863         return
 864
 865     def _as_decimal(self):
 866         return decimal.Decimal('-Infinity')
 867
 868 class NaN(Singleton, Rational):
 869
 870     p = 0
 871     q = 0
 872
 873     is_commutative = True
 874     is_real = None
 875     is_comparable = None
 876     is_bounded = None
 877     #is_unbounded = False
 878
 879     def tostr(self, level=0):
 880         return 'nan'
 881
 882     def _as_decimal(self):
 883         return decimal.Decimal('NaN')
 884
 885     def _eval_power(b, e):
 886         if isinstance(e, Basic.Zero):
 887             return S.One
 888         return b
 889
 890 class ComplexInfinity(Singleton, Atom, NoRelMeths, ArithMeths):
 891
 892     is_commutative = True
 893     is_comparable = None
 894     is_bounded = False
 895     is_real = None
 896
 897     def tostr(self, level=0):
 898         return 'zoo'
 899
 900     def _eval_power(b, e):
 901         if isinstance(e, Basic.ComplexInfinity):
 902             return S.NaN
 903
 904         if isinstance(e, Basic.Number):
 905             if isinstance(e, Basic.Zero):
 906                 return S.NaN
 907             else:
 908                 if e.is_positive:
 909                     return S.ComplexInfinity
 910                 else:
 911                     return S.Zero
 912
 913 class NumberSymbol(Singleton, Atom, RelMeths, ArithMeths):
 914
 915     is_commutative = True
 916     is_comparable = True
 917     is_bounded = True
 918     is_finite = True
 919
 920     def approximation(self, number_cls):
 921         """ Return an interval with number_cls endpoints
 922         that contains the value of NumberSymbol.
 923         If not implemented, then return None.
 924         """
 925
 926     def _eval_derivative(self, s):
 927         return S.Zero
 928     def __eq__(self, other):
 929         other = Basic.sympify(other)
 930         if self is other: return True
 931         if isinstance(other, Number) and self.is_irrational: return False
 932         return RelMeths.__eq__(self, other)
 933     def __ne__(self, other):
 934         other = Basic.sympify(other)
 935         if self is other: return False
 936         if isinstance(other, Number) and self.is_irrational: return True
 937         return RelMeths.__ne__(self, other)
 938     def __lt__(self, other):
 939         other = Basic.sympify(other)
 940         if self is other: return False
 941         if isinstance(other, Number):
 942             approx = self.approximation_interval(other.__class__)
 943             if approx is not None:
 944                 l,u = approx
 945                 if other < l: return False
 946                 if other > u: return True
 947             return self.evalf()<other
 948         if other.is_comparable:
 949             other = other.evalf()
 950             return self.evalf()<other
 951         return RelMeths.__lt__(self, other)
 952     def __le__(self, other):
 953         other = Basic.sympify(other)
 954         if self is other: return True
 955         if other.is_comparable: other = other.evalf()
 956         if isinstance(other, Number):
 957             return self.evalf()<=other
 958         return RelMeths.__le__(self, other)
 959     def __gt__(self, other):
 960         return (-self) < (-other)
 961     def __ge__(self, other):
 962         return (-self) <= (-other)
 963
 964
 965 class Exp1(NumberSymbol):
 966
 967     is_real = True
 968     is_positive = True
 969     is_negative = False # XXX Forces is_negative/is_nonnegative
 970     is_irrational = True
 971
 972     def tostr(self, level=0):
 973         return 'E'
 974
 975     def _eval_evalf(self):
 976         return Real(decimal_math.e())
 977
 978     def approximation_interval(self, number_cls):
 979         if issubclass(number_cls,Integer):
 980             return (Integer(2),Integer(3))
 981         elif issubclass(number_cls,Rational):
 982             pass
 983
 984     def _eval_power(self, exp):
 985         return Basic.exp(exp)
 986
 987 class Pi(NumberSymbol):
 988
 989     is_real = True
 990     is_positive = True
 991     is_negative = False
 992     is_irrational = True
 993
 994     def _eval_evalf(self):
 995         return Real(decimal_math.pi())
 996
 997     def approximation_interval(self, number_cls):
 998         if issubclass(number_cls, Integer):
 999             return (Integer(3), Integer(4))
1000         elif issubclass(number_cls, Rational):
1001             return (Rational(223,71), Rational(22,7))
1002
1003     def tostr(self, level=0):
1004         return 'pi'
1005
1006 class GoldenRatio(NumberSymbol):
1007
1008     is_real = True
1009     is_positive = True
1010     is_negative = False
1011     is_irrational = True
1012
1013     def _eval_evalf(self):
1014         return Real(decimal_math.golden_ratio())
1015
1016     def _eval_expand_func(self, *args):
1017         return S.Half + S.Half*S.Sqrt(5)
1018
1019     def approximation_interval(self, number_cls):
1020         if issubclass(number_cls, Integer):
1021             return (S.One, Rational(2))
1022         elif issubclass(number_cls, Rational):
1023             pass
1024
1025     def tostr(self, level=0):
1026         return 'GoldenRatio'
1027
1028 class EulerGamma(NumberSymbol):
1029
1030     is_real = True
1031     is_positive = True
1032     is_negative = False
1033     is_irrational = None
1034
1035     def _eval_evalf(self):
1036         return
1037
1038     def approximation_interval(self, number_cls):
1039         if issubclass(number_cls, Integer):
1040             return (S.Zero, S.One)
1041         elif issubclass(number_cls, Rational):
1042             return (S.Half, Rational(3, 5))
1043
1044     def tostr(self, level=0):
1045         return 'EulerGamma'
1046
1047 class Catalan(NumberSymbol):
1048
1049     is_real = True
1050     is_positive = True
1051     is_negative = False
1052     is_irrational = None
1053
1054     def _eval_evalf(self):
1055         return
1056
1057     def approximation_interval(self, number_cls):
1058         if issubclass(number_cls, Integer):
1059             return (S.Zero, S.One)
1060         elif issubclass(number_cls, Rational):
1061             return (Rational(9, 10), S.One)
1062
1063     def tostr(self, level=0):
1064         return 'Catalan'
1065
1066 class ImaginaryUnit(Singleton, Atom, RelMeths, ArithMeths):
1067
1068     is_commutative = True
1069     is_imaginary = True
1070     is_bounded = True
1071     is_finite = True
1072
1073     def tostr(self, level=0):
1074         return 'I'
1075
1076     def _eval_conjugate(self):
1077         return -S.ImaginaryUnit
1078
1079     def _eval_derivative(self, s):
1080         return S.Zero
1081
1082     def _eval_power(b, e):
1083         """
1084         b is I = sqrt(-1)
1085         e is symbolic object but not equal to 0, 1
1086
1087         I ** r -> (-1)**(r/2) -> exp(r/2 * Pi * I) -> sin(Pi*r/2) + cos(Pi*r/2) * I, r is decimal
1088         I ** 0 mod 4 -> 1
1089         I ** 1 mod 4 -> I
1090         I ** 2 mod 4 -> -1
1091         I ** 3 mod 4 -> -I
1092         """
1093
1094
1095         if isinstance(e, Number):
1096             #if isinstance(e, Decimal):
1097             #    a = decimal_math.pi() * exponent.num / 2
1098             #    return Decimal(decimal_math.sin(a) + decimal_math.cos(a) * ImaginaryUnit())
1099             if isinstance(e, Integer):
1100                 e = e.p % 4
1101                 if e==0: return S.One
1102                 if e==1: return S.ImaginaryUnit
1103                 if e==2: return -S.One
1104                 return -S.ImaginaryUnit
1105             return (-S.One) ** (e * S.Half)
1106         return
1107
1108     def as_base_exp(self):
1109         return -S.One, S.Half
1110
1111 Basic.singleton['E'] = Exp1
1112 Basic.singleton['pi'] = Pi
1113 Basic.singleton['I'] = ImaginaryUnit
1114 Basic.singleton['oo'] = Infinity
1115 Basic.singleton['nan'] = NaN
1116
1117 Basic.singleton['zoo'] = ComplexInfinity
1118
1119 Basic.singleton['GoldenRatio'] = GoldenRatio
1120 Basic.singleton['EulerGamma'] = EulerGamma
1121 Basic.singleton['Catalan'] = Catalan