sympy/core/power.py

   1
   2 from basic import Basic, S, cache_it, cache_it_immutable
   3 from methods import ArithMeths, RelMeths, NoRelMeths
   4
   5
   6 def integer_nthroot(y, n):
   7     """
   8     Usage
   9     =====
  10         Return a tuple containing x = floor(y**(1/n))
  11         and a boolean indicating whether the result is exact (that is,
  12         whether x**n == y).
  13
  14     Examples
  15     ========
  16
  17         >>> integer_nthroot(16,2)
  18         (4, True)
  19         >>> integer_nthroot(26,2)
  20         (5, False)
  21         >>> integer_nthroot(1234567**7, 7)
  22         (1234567L, True)
  23         >>> integer_nthroot(1234567**7+1, 7)
  24         (1234567L, False)
  25     """
  26
  27     y = int(y); n = int(n)
  28     if y < 0: raise ValueError, "y must not be negative"
  29     if n < 1: raise ValueError, "n must be positive"
  30     if y in (0, 1): return y, True
  31     if n == 1: return y, True
  32
  33     # Search with Newton's method, starting from floating-point
  34     # approximation. Care must be taken to avoid overflow.
  35     from math import log as _log
  36     guess = 2**int(_log(y, 2)/n)
  37     xprev, x = -1, guess
  38     while abs(x - xprev) > 1:
  39         t = x**(n-1)
  40         xprev, x = x, x - (t*x-y)//(n*t)
  41     # Compensate
  42     while x**n > y:
  43         x -= 1
  44     return x, x**n == y
  45
  46
  47 class Pow(Basic, ArithMeths, RelMeths):
  48
  49     precedence = Basic.Pow_precedence
  50
  51     @cache_it_immutable
  52     def __new__(cls, a, b, **assumptions):
  53         a = Basic.sympify(a)
  54         b = Basic.sympify(b)
  55         if isinstance(b, Basic.Zero):
  56             return S.One
  57         if isinstance(b, Basic.One):
  58             return a
  59         obj = a._eval_power(b)
  60         if obj is None:
  61             obj = Basic.__new__(cls, a, b, **assumptions)
  62         return obj
  63
  64     @property
  65     def base(self):
  66         return self._args[0]
  67
  68     @property
  69     def exp(self):
  70         return self._args[1]
  71
  72     def _eval_power(self, other):
  73         if isinstance(other, Basic.Number):
  74             if self.base.is_real:
  75                 if isinstance(self.exp, Basic.Number):
  76                     # (a ** 2) ** 3 -> a ** (2 * 3)
  77                     return Pow(self.base, self.exp * other)
  78             if isinstance(other, Basic.Integer):
  79                 # (a ** b) ** 3 -> a ** (3 * b)
  80                 return Pow(self.base, self.exp * other)
  81         elif isinstance(other, (Basic.Add, Basic.Mul)):
  82             # (a**b)**c = a**(b*c)
  83             return Pow(self.base, self.exp * other)
  84
  85         if other.atoms(Basic.Wild):
  86             return Pow(self.base, self.exp * other)
  87         return
  88
  89     def _eval_is_commutative(self):
  90         c1 = self.base.is_commutative
  91         if c1 is None: return
  92         c2 = self.base.is_commutative
  93         if c2 is None: return
  94         return c1 and c2
  95
  96     def _eval_is_comparable(self):
  97         c1 = self.base.is_comparable
  98         if c1 is None: return
  99         c2 = self.base.is_comparable
 100         if c2 is None: return
 101         return c1 and c2
 102
 103     def _eval_is_even(self):
 104         if self.exp.is_integer and self.exp.is_positive:
 105             if self.base.is_even:
 106                 return True
 107             if self.exp.is_integer:
 108                 return False
 109
 110     def _eval_is_positive(self):
 111         if self.base.is_positive:
 112             if self.exp.is_real:
 113                 return True
 114         elif self.base.is_negative:
 115             if self.exp.is_even:
 116                 return True
 117             if self.exp.is_odd:
 118                 return False
 119         elif self.base.is_nonpositive:
 120             if self.exp.is_odd:
 121                 return False
 122
 123     def _eval_is_negative(self):
 124         if self.base.is_negative:
 125             if self.exp.is_odd:
 126                 return True
 127             if self.exp.is_even:
 128                 return False
 129         elif self.base.is_positive:
 130             if self.exp.is_real:
 131                 return False
 132         elif self.base.is_nonnegative:
 133             if self.exp.is_real:
 134                 return False
 135         elif self.base.is_nonpositive:
 136             if self.exp.is_even:
 137                 return False
 138         elif self.base.is_real:
 139             if self.exp.is_even:
 140                 return False
 141
 142     def _eval_is_integer(self):
 143         c1 = self.base.is_integer
 144         if c1 is None: return
 145         c2 = self.exp.is_integer
 146         if c2 is None: return
 147         if c1 and c2:
 148             if self.exp.is_nonnegative or self.exp.is_positive:
 149                 return True
 150             if self.exp.is_negative:
 151                 return False
 152
 153     def _eval_is_real(self):
 154         c1 = self.base.is_real
 155         if c1 is None: return
 156         c2 = self.exp.is_real
 157         if c2 is None: return
 158         if c1 and c2:
 159             if self.base.is_positive:
 160                 return True
 161             if self.base.is_negative:
 162                 if self.exp.is_integer:
 163                     return True
 164
 165     def _eval_is_odd(self):
 166         if not (self.base.is_integer and self.exp.is_nonnegative): return
 167         return self.base.is_odd
 168
 169     def _eval_is_bounded(self):
 170         if self.exp.is_negative:
 171             if self.base.is_infinitesimal:
 172                 return False
 173             if self.base.is_unbounded:
 174                 return True
 175         c1 = self.base.is_bounded
 176         if c1 is None: return
 177         c2 = self.exp.is_bounded
 178         if c2 is None: return
 179         if c1 and c2:
 180             if self.exp.is_nonnegative:
 181                 return True
 182
 183     def tostr(self, level=0):
 184         precedence = self.precedence
 185         b = self.base.tostr(precedence)
 186         if isinstance(self.exp, Basic.NegativeOne):
 187             r = '1/%s' % (b)
 188         else:
 189             r = '%s**%s' % (b,self.exp.tostr(precedence))
 190         if precedence <= level:
 191             return '(%s)' % (r)
 192         return r
 193
 194     def _eval_subs(self, old, new):
 195         if self==old: return new
 196         if isinstance(old, self.__class__) and self.base==old.base:
 197             coeff1,terms1 = self.exp.as_coeff_terms()
 198             coeff2,terms2 = old.exp.as_coeff_terms()
 199             if terms1==terms2: return new ** (coeff1/coeff2) # (x**(2*y)).subs(x**(3*y),z) -> z**(2/3*y)
 200         if isinstance(old, Basic.ApplyExp):
 201             coeff1,terms1 = old[0].as_coeff_terms()
 202             coeff2,terms2 = (self.exp * S.Log(self.base)).as_coeff_terms()
 203             if terms1==terms2: return new ** (coeff1/coeff2) # (x**(2*y)).subs(exp(3*y*log(x)),z) -> z**(2/3*y)
 204         return self.base.subs(old, new) ** self.exp.subs(old, new)
 205
 206     def as_powers_dict(self):
 207         return { self.base : self.exp }
 208
 209     def as_base_exp(self):
 210         if isinstance(self.base, Basic.Rational) and self.base.p==1:
 211             return 1/self.base, -self.exp
 212         return self.base, self.exp
 213
 214     def _eval_conjugate(self):
 215         return S.Conjugate(self.base)**self.exp
 216
 217     def _eval_expand_complex(self, *args):
 218         if isinstance(self.exp, Basic.Integer):
 219             re, im = self.base.as_real_imag()
 220             base = re + S.ImaginaryUnit * im
 221             return (base**self.exp)._eval_expand_basic(*args)
 222         elif isinstance(self.exp, Basic.Rational):
 223             # NOTE: This is not totally correct since for x**(p/q) with
 224             #       x being imaginary there are actually q roots, but
 225             #       only a single one is returned from here.
 226             re, im = self.base.as_real_imag()
 227
 228             r = S.Sqrt(re**2 + im**2)
 229             t = Basic.atan(im / re)
 230
 231             if im == 0 and re == -1:
 232                 t = S.Pi
 233
 234             rp, tp = r**self.exp, t*self.exp
 235
 236             return rp*Basic.cos(tp) + rp*Basic.sin(tp)*S.ImaginaryUnit
 237         else:
 238             return S.Re(self) + S.ImaginaryUnit*S.Im(self)
 239
 240     def _eval_expand_basic(self, *args):
 241         """
 242         (a*b)**n -> a**n * b**n
 243         (a+b+..) ** n -> a**n + n*a**(n-1)*b + .., n is positive integer
 244         """
 245         base = self.base.expand()
 246         exponent = self.exp.expand()
 247         result = base ** exponent
 248         if isinstance(result, Pow):
 249             base = result.base
 250             exponent = result.exp
 251         else:
 252             return result
 253         if isinstance(exponent, Basic.Integer):
 254             if isinstance(base, Basic.Mul):
 255                 return Basic.Mul(*[t**exponent for t in base])
 256             if exponent.is_positive and isinstance(base, Basic.Add):
 257                 m = int(exponent)
 258                 if base.is_commutative:
 259                     p = []
 260                     order_terms = []
 261                     for o in base:
 262                         if isinstance(o, Basic.Order):
 263                             order_terms.append(o)
 264                         else:
 265                             p.append(o)
 266                     if order_terms:
 267                         # (f(x) + O(x^n))^m -> f(x)^m + m*f(x)^{m-1} *O(x^n)
 268                         f = Basic.Add(*p)
 269                         fm1 = (f**(m-1)).expand()
 270                         return (f*fm1).expand() + m*fm1*Basic.Add(*order_terms)
 271                 ## Consider polynomial
 272                 ##   P(x) = sum_{i=0}^n p_i x^k
 273                 ## and its m-th exponent
 274                 ##   P(x)^m = sum_{k=0}^{m n} a(m,k) x^k
 275                 ## The coefficients a(m,k) can be computed using the
 276                 ## J.C.P. Miller Pure Recurrence [see D.E.Knuth,
 277                 ## Seminumerical Algorithms, The art of Computer
 278                 ## Programming v.2, Addison Wesley, Reading, 1981;]:
 279                 ##  a(m,k) = 1/(k p_0) sum_{i=1}^n p_i ((m+1)i-k) a(m,k-i),
 280                 ## where a(m,0) = p_0^m.
 281                     n = len(p)-1
 282                     cache = {0: p[0] ** m}
 283                     p0 = [t/p[0] for t in p]
 284                     l = [cache[0]]
 285                     Mul = Basic.Mul
 286                     Rational = Basic.Rational
 287                     for k in xrange(1, m * n + 1):
 288                         a = []
 289                         for i in xrange(1,n+1):
 290                             if i<=k:
 291                                 a.append(Mul(Rational((m+1)*i-k,k), p0[i], cache[k-i]).expand())
 292                         a = Basic.Add(*a)
 293                         cache[k] = a
 294                         l.append(a)
 295                     return Basic.Add(*l)
 296                 else:
 297                     if m==2:
 298                         p = base[:]
 299                         return Basic.Add(*[t1*t2 for t1 in p for t2 in p])
 300                     return Basic.Mul(base, Pow(base, m-1).expand()).expand()
 301         elif isinstance(exponent, Basic.Add) and isinstance(base, Basic.Number):
 302             # n**(a+b) --> n**a * n**b with n and a Numbers
 303             exp = 0
 304             coeff = 1
 305             for term in exponent:
 306                 if isinstance(term, Basic.Number):
 307                     coeff *= base**term
 308                 else:
 309                     exp += term
 310             result = coeff * base**exp
 311
 312         return result
 313
 314     def _eval_derivative(self, s):
 315         dbase = self.base.diff(s)
 316         dexp = self.exp.diff(s)
 317         return self * (dexp * S.Log(self.base) + dbase * self.exp/self.base)
 318
 319     _eval_evalf = Basic._seq_eval_evalf
 320
 321     def _calc_splitter(self, d):
 322         if d.has_key(self):
 323             return d[self]
 324         base = self.base._calc_splitter(d)
 325         exp = self.exp._calc_splitter(d)
 326         if isinstance(exp, Basic.Integer):
 327             if abs(exp.p)>2:
 328                 n = exp.p//2
 329                 r = exp.p - n
 330                 if n!=r:
 331                     p1 = (base ** n)._calc_splitter(d)
 332                     p2 = (base ** r)._calc_splitter(d)
 333                     r = p1*p2
 334                 else:
 335                     r = (base ** n)._calc_splitter(d) ** 2
 336             elif exp.p==-2:
 337                 r = (1/base)._calc_splitter(d) ** 2
 338             else:
 339                 r = base ** exp
 340         else:
 341             r = base ** exp
 342         if d.has_key(r):
 343             return d[r]
 344         s = d[r] = Basic.Temporary()
 345         return s
 346
 347     @cache_it_immutable
 348     def count_ops(self, symbolic=True):
 349         if symbolic:
 350             return Basic.Add(*[t.count_ops(symbolic) for t in self[:]]) + Basic.Symbol('POW')
 351         return Basic.Add(*[t.count_ops(symbolic) for t in self[:]]) + 1
 352
 353     def _eval_integral(self, s):
 354         if not self.exp.has(s):
 355             if self.base==s:
 356                 n = self.exp+1
 357                 return self.base ** n/n
 358             y = Basic.Symbol('y',dummy=True)
 359             x,ix = self.base.solve4linearsymbol(y,symbols=set([s]))
 360             if isinstance(x, Basic.Symbol):
 361                 dx = 1/self.base.diff(x)
 362                 if not dx.has(s):
 363                     return (y**self.exp*dx).integral(y).subs(y, self.base)
 364
 365     def _eval_defined_integral(self, s, a, b):
 366         if not self.exp.has(s):
 367             if self.base==s:
 368                 n = self.exp+1
 369                 return (b**n-a**n)/n
 370             x,ix = self.base.solve4linearsymbol(s)
 371             if isinstance(x, Basic.Symbol):
 372                 dx = ix.diff(x)
 373                 if isinstance(dx, Basic.Number):
 374                     y = Basic.Symbol('y',dummy=True)
 375                     return (y**self.exp*dx).integral(y==[self.base.subs(s,a), self.base.subs(s,b)])
 376
 377     def _eval_is_polynomial(self, syms):
 378         if self.exp.has(*syms):
 379             return False
 380
 381         if self.base.has(*syms):
 382             # it would be nice to have is_nni working
 383             return self.base._eval_is_polynomial(syms) and \
 384                    self.exp.is_nonnegative and \
 385                    self.exp.is_integer
 386         else:
 387             return True
 388
 389     def as_numer_denom(self):
 390         base, exp = self.as_base_exp()
 391         c,t = exp.as_coeff_terms()
 392         n,d = base.as_numer_denom()
 393         if c.is_negative:
 394             exp = -exp
 395             n,d = d,n
 396         return n ** exp, d ** exp
 397
 398     def matches(pattern, expr, repl_dict={}, evaluate=False):
 399         Basic.matches.__doc__
 400         if evaluate:
 401             pat = pattern
 402             for old,new in repl_dict.items():
 403                 pat = pat.subs(old, new)
 404             if pat!=pattern:
 405                 return pat.matches(expr, repl_dict)
 406
 407         expr = Basic.sympify(expr)
 408         b, e = expr.as_base_exp()
 409
 410         # special case, pattern = 1 and expr.exp can match to 0
 411         if isinstance(expr, Basic.One):
 412             d = repl_dict.copy()
 413             d = pattern.exp.matches(Basic.Integer(0), d, evaluate=False)
 414             if d is not None:
 415                 return d
 416
 417         d = repl_dict.copy()
 418         d = pattern.base.matches(b, d, evaluate=False)
 419         if d is None:
 420             return None
 421
 422         d = pattern.exp.matches(e, d, evaluate=True)
 423         if d is None:
 424             return Basic.matches(pattern, expr, repl_dict, evaluate)
 425         return d
 426
 427     def _eval_oseries(self, order):
 428         """
 429         f**g + O(h) == (f+O(k))**g -> .. -> f**g + O(f**(g-1)*k), hence O(k)==O(h*f**(1-g)).
 430         If f->0 as x->0 then
 431         """
 432         x = order.symbols[0]
 433         e = self.exp
 434         b = self.base
 435         ln = S.Log
 436         exp = S.Exp
 437         if e.has(x):
 438             return exp(e * ln(b)).oseries(order)
 439         if b==x: return self
 440         b0 = b.limit(x,0)
 441         if isinstance(b0, Basic.Zero) or b0.is_unbounded:
 442             lt = b.as_leading_term(x)
 443             o = order * lt**(1-e)
 444             bs = b.oseries(o)
 445             if isinstance(bs, Basic.Add):
 446                 # bs -> lt + rest -> lt * (1 + (bs/lt - 1))
 447                 return lt**e * ((bs/lt).expand()**e).oseries(order * lt**(-e))
 448             return bs**e
 449         o = order * (b0**-e)
 450         # b -> b0 + (b-b0) -> b0 * (1 + (b/b0-1))
 451         z = (b/b0-1)
 452         r = self._compute_oseries(z, o, self.taylor_term, lambda z: 1+z) * b0**e
 453         return r
 454
 455     def _eval_as_leading_term(self, x):
 456         if not self.exp.has(x):
 457             return self.base.as_leading_term(x) ** self.exp
 458         return S.Exp(self.exp * S.Log(self.base)).as_leading_term(x)
 459
 460     @cache_it_immutable
 461     def taylor_term(self, n, x, *previous_terms): # of (1+x)**e
 462         if n<0: return S.Zero
 463         x = Basic.sympify(x)
 464         return S.Binomial(self.exp, n) * x**n