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.exp):
 201             coeff1,terms1 = old[0].as_coeff_terms()
 202             coeff2,terms2 = (self.exp * Basic.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         from sympy.functions.elementary.complexes import conjugate as c
 216         return c(self.base)**self.exp
 217
 218     def _eval_expand_complex(self, *args):
 219         if isinstance(self.exp, Basic.Integer):
 220             re, im = self.base.as_real_imag()
 221             base = re + S.ImaginaryUnit * im
 222             return (base**self.exp)._eval_expand_basic(*args)
 223         elif isinstance(self.exp, Basic.Rational):
 224             # NOTE: This is not totally correct since for x**(p/q) with
 225             #       x being imaginary there are actually q roots, but
 226             #       only a single one is returned from here.
 227             re, im = self.base.as_real_imag()
 228
 229             r = S.Sqrt(re**2 + im**2)
 230             t = Basic.atan(im / re)
 231
 232             if im == 0 and re == -1:
 233                 t = S.Pi
 234
 235             rp, tp = r**self.exp, t*self.exp
 236
 237             return rp*Basic.cos(tp) + rp*Basic.sin(tp)*S.ImaginaryUnit
 238         else:
 239             return S.Re(self) + S.ImaginaryUnit*S.Im(self)
 240
 241     def _eval_expand_basic(self, *args):
 242         """
 243         (a*b)**n -> a**n * b**n
 244         (a+b+..) ** n -> a**n + n*a**(n-1)*b + .., n is positive integer
 245         """
 246         base = self.base.expand()
 247         exponent = self.exp.expand()
 248         result = base ** exponent
 249         if isinstance(result, Pow):
 250             base = result.base
 251             exponent = result.exp
 252         else:
 253             return result
 254         if isinstance(exponent, Basic.Integer):
 255             if isinstance(base, Basic.Mul):
 256                 return Basic.Mul(*[t**exponent for t in base])
 257             if exponent.is_positive and isinstance(base, Basic.Add):
 258                 m = int(exponent)
 259                 if base.is_commutative:
 260                     p = []
 261                     order_terms = []
 262                     for o in base:
 263                         if isinstance(o, Basic.Order):
 264                             order_terms.append(o)
 265                         else:
 266                             p.append(o)
 267                     if order_terms:
 268                         # (f(x) + O(x^n))^m -> f(x)^m + m*f(x)^{m-1} *O(x^n)
 269                         f = Basic.Add(*p)
 270                         fm1 = (f**(m-1)).expand()
 271                         return (f*fm1).expand() + m*fm1*Basic.Add(*order_terms)
 272                 ## Consider polynomial
 273                 ##   P(x) = sum_{i=0}^n p_i x^k
 274                 ## and its m-th exponent
 275                 ##   P(x)^m = sum_{k=0}^{m n} a(m,k) x^k
 276                 ## The coefficients a(m,k) can be computed using the
 277                 ## J.C.P. Miller Pure Recurrence [see D.E.Knuth,
 278                 ## Seminumerical Algorithms, The art of Computer
 279                 ## Programming v.2, Addison Wesley, Reading, 1981;]:
 280                 ##  a(m,k) = 1/(k p_0) sum_{i=1}^n p_i ((m+1)i-k) a(m,k-i),
 281                 ## where a(m,0) = p_0^m.
 282                     n = len(p)-1
 283                     cache = {0: p[0] ** m}
 284                     p0 = [t/p[0] for t in p]
 285                     l = [cache[0]]
 286                     Mul = Basic.Mul
 287                     Rational = Basic.Rational
 288                     for k in xrange(1, m * n + 1):
 289                         a = []
 290                         for i in xrange(1,n+1):
 291                             if i<=k:
 292                                 a.append(Mul(Rational((m+1)*i-k,k), p0[i], cache[k-i]).expand())
 293                         a = Basic.Add(*a)
 294                         cache[k] = a
 295                         l.append(a)
 296                     return Basic.Add(*l)
 297                 else:
 298                     if m==2:
 299                         p = base[:]
 300                         return Basic.Add(*[t1*t2 for t1 in p for t2 in p])
 301                     return Basic.Mul(base, Pow(base, m-1).expand()).expand()
 302         elif isinstance(exponent, Basic.Add) and isinstance(base, Basic.Number):
 303             # n**(a+b) --> n**a * n**b with n and a Numbers
 304             exp = 0
 305             coeff = 1
 306             for term in exponent:
 307                 if isinstance(term, Basic.Number):
 308                     coeff *= base**term
 309                 else:
 310                     exp += term
 311             result = coeff * base**exp
 312
 313         return result
 314
 315     def _eval_derivative(self, s):
 316         dbase = self.base.diff(s)
 317         dexp = self.exp.diff(s)
 318         return self * (dexp * Basic.log(self.base) + dbase * self.exp/self.base)
 319
 320     _eval_evalf = Basic._seq_eval_evalf
 321
 322     def _calc_splitter(self, d):
 323         if d.has_key(self):
 324             return d[self]
 325         base = self.base._calc_splitter(d)
 326         exp = self.exp._calc_splitter(d)
 327         if isinstance(exp, Basic.Integer):
 328             if abs(exp.p)>2:
 329                 n = exp.p//2
 330                 r = exp.p - n
 331                 if n!=r:
 332                     p1 = (base ** n)._calc_splitter(d)
 333                     p2 = (base ** r)._calc_splitter(d)
 334                     r = p1*p2
 335                 else:
 336                     r = (base ** n)._calc_splitter(d) ** 2
 337             elif exp.p==-2:
 338                 r = (1/base)._calc_splitter(d) ** 2
 339             else:
 340                 r = base ** exp
 341         else:
 342             r = base ** exp
 343         if d.has_key(r):
 344             return d[r]
 345         s = d[r] = Basic.Temporary()
 346         return s
 347
 348     @cache_it_immutable
 349     def count_ops(self, symbolic=True):
 350         if symbolic:
 351             return Basic.Add(*[t.count_ops(symbolic) for t in self[:]]) + Basic.Symbol('POW')
 352         return Basic.Add(*[t.count_ops(symbolic) for t in self[:]]) + 1
 353
 354     def _eval_integral(self, s):
 355         if not self.exp.has(s):
 356             if self.base==s:
 357                 n = self.exp+1
 358                 return self.base ** n/n
 359             y = Basic.Symbol('y',dummy=True)
 360             x,ix = self.base.solve4linearsymbol(y,symbols=set([s]))
 361             if isinstance(x, Basic.Symbol):
 362                 dx = 1/self.base.diff(x)
 363                 if not dx.has(s):
 364                     return (y**self.exp*dx).integral(y).subs(y, self.base)
 365
 366     def _eval_defined_integral(self, s, a, b):
 367         if not self.exp.has(s):
 368             if self.base==s:
 369                 n = self.exp+1
 370                 return (b**n-a**n)/n
 371             x,ix = self.base.solve4linearsymbol(s)
 372             if isinstance(x, Basic.Symbol):
 373                 dx = ix.diff(x)
 374                 if isinstance(dx, Basic.Number):
 375                     y = Basic.Symbol('y',dummy=True)
 376                     return (y**self.exp*dx).integral(y==[self.base.subs(s,a), self.base.subs(s,b)])
 377
 378     def _eval_is_polynomial(self, syms):
 379         if self.exp.has(*syms):
 380             return False
 381
 382         if self.base.has(*syms):
 383             # it would be nice to have is_nni working
 384             return self.base._eval_is_polynomial(syms) and \
 385                    self.exp.is_nonnegative and \
 386                    self.exp.is_integer
 387         else:
 388             return True
 389
 390     def as_numer_denom(self):
 391         base, exp = self.as_base_exp()
 392         c,t = exp.as_coeff_terms()
 393         n,d = base.as_numer_denom()
 394         if c.is_negative:
 395             exp = -exp
 396             n,d = d,n
 397         return n ** exp, d ** exp
 398
 399     def matches(pattern, expr, repl_dict={}, evaluate=False):
 400         Basic.matches.__doc__
 401         if evaluate:
 402             pat = pattern
 403             for old,new in repl_dict.items():
 404                 pat = pat.subs(old, new)
 405             if pat!=pattern:
 406                 return pat.matches(expr, repl_dict)
 407
 408         expr = Basic.sympify(expr)
 409         b, e = expr.as_base_exp()
 410
 411         # special case, pattern = 1 and expr.exp can match to 0
 412         if isinstance(expr, Basic.One):
 413             d = repl_dict.copy()
 414             d = pattern.exp.matches(Basic.Integer(0), d, evaluate=False)
 415             if d is not None:
 416                 return d
 417
 418         d = repl_dict.copy()
 419         d = pattern.base.matches(b, d, evaluate=False)
 420         if d is None:
 421             return None
 422
 423         d = pattern.exp.matches(e, d, evaluate=True)
 424         if d is None:
 425             return Basic.matches(pattern, expr, repl_dict, evaluate)
 426         return d
 427
 428     def _eval_oseries(self, order):
 429         """
 430         f**g + O(h) == (f+O(k))**g -> .. -> f**g + O(f**(g-1)*k), hence O(k)==O(h*f**(1-g)).
 431         If f->0 as x->0 then
 432         """
 433         x = order.symbols[0]
 434         e = self.exp
 435         b = self.base
 436         ln = Basic.log
 437         exp = Basic.exp
 438         if e.has(x):
 439             return exp(e * ln(b)).oseries(order)
 440         if b==x: return self
 441         b0 = b.limit(x,0)
 442         if isinstance(b0, Basic.Zero) or b0.is_unbounded:
 443             lt = b.as_leading_term(x)
 444             o = order * lt**(1-e)
 445             bs = b.oseries(o)
 446             if isinstance(bs, Basic.Add):
 447                 # bs -> lt + rest -> lt * (1 + (bs/lt - 1))
 448                 return lt**e * ((bs/lt).expand()**e).oseries(order * lt**(-e))
 449             return bs**e
 450         o = order * (b0**-e)
 451         # b -> b0 + (b-b0) -> b0 * (1 + (b/b0-1))
 452         z = (b/b0-1)
 453         r = self._compute_oseries(z, o, self.taylor_term, lambda z: 1+z) * b0**e
 454         return r
 455
 456     def _eval_as_leading_term(self, x):
 457         if not self.exp.has(x):
 458             return self.base.as_leading_term(x) ** self.exp
 459         return Basic.exp(self.exp * Basic.log(self.base)).as_leading_term(x)
 460
 461     @cache_it_immutable
 462     def taylor_term(self, n, x, *previous_terms): # of (1+x)**e
 463         if n<0: return S.Zero
 464         x = Basic.sympify(x)
 465         return S.Binomial(self.exp, n) * x**n