1import math 2 3from mpmath.libmp import sqrtrem as mpmath_sqrtrem 4 5from ..logic import true 6from ..utilities import sift 7from .add import Add 8from .cache import cacheit 9from .compatibility import as_int 10from .evalf import PrecisionExhausted 11from .evaluate import global_evaluate 12from .expr import Expr 13from .function import (_coeff_isneg, expand_complex, expand_mul, 14 expand_multinomial) 15from .logic import fuzzy_or 16from .mul import Mul, _keep_coeff 17from .numbers import E, I, Integer, Rational, nan, oo, pi, zoo 18from .singleton import S 19from .symbol import Dummy, symbols 20from .sympify import sympify 21 22 23def integer_nthroot(y, n): 24 """ 25 Return a tuple containing x = floor(y**(1/n)) 26 and a boolean indicating whether the result is exact (that is, 27 whether x**n == y). 28 29 >>> integer_nthroot(16, 2) 30 (4, True) 31 >>> integer_nthroot(26, 2) 32 (5, False) 33 34 """ 35 y, n = int(y), int(n) 36 if y < 0: 37 raise ValueError('y must be nonnegative') 38 if n < 1: 39 raise ValueError('n must be positive') 40 if y in (0, 1): 41 return y, True 42 if n == 1: 43 return y, True 44 if n == 2: 45 x, rem = mpmath_sqrtrem(y) 46 return int(x), not rem 47 if n > y: 48 return 1, False 49 # Get initial estimate for Newton's method. Care must be taken to 50 # avoid overflow 51 try: 52 guess = int(y**(1./n) + 0.5) 53 except OverflowError: 54 exp = math.log(y, 2)/n 55 if exp > 53: 56 shift = int(exp - 53) 57 guess = int(2.0**(exp - shift) + 1) << shift 58 else: 59 guess = int(2.0**exp) 60 if guess > 2**50: 61 # Newton iteration 62 xprev, x = -1, guess 63 while 1: 64 t = x**(n - 1) 65 xprev, x = x, ((n - 1)*x + y//t)//n 66 if abs(x - xprev) < 2: 67 break 68 else: 69 x = guess 70 # Compensate 71 t = x**n 72 while t < y: 73 x += 1 74 t = x**n 75 while t > y: 76 x -= 1 77 t = x**n 78 return x, t == y 79 80 81class Pow(Expr): 82 """ 83 Defines the expression x**y as "x raised to a power y". 84 85 For complex numbers `x` and `y`, ``Pow`` gives the principal 86 value of `exp(y*log(x))`. 87 88 Singleton definitions involving (0, 1, -1, oo, -oo, I, -I): 89 90 +--------------+---------+-----------------------------------------------+ 91 | expr | value | reason | 92 +==============+=========+===============================================+ 93 | z**0 | 1 | Although arguments over 0**0 exist, see [2]. | 94 +--------------+---------+-----------------------------------------------+ 95 | z**1 | z | | 96 +--------------+---------+-----------------------------------------------+ 97 | (-oo)**(-1) | 0 | | 98 +--------------+---------+-----------------------------------------------+ 99 | (-1)**-1 | -1 | | 100 +--------------+---------+-----------------------------------------------+ 101 | 0**-1 | zoo | This is not strictly true, as 0**-1 may be | 102 | | | undefined, but is convenient in some contexts | 103 | | | where the base is assumed to be positive. | 104 +--------------+---------+-----------------------------------------------+ 105 | 1**-1 | 1 | | 106 +--------------+---------+-----------------------------------------------+ 107 | oo**-1 | 0 | | 108 +--------------+---------+-----------------------------------------------+ 109 | 0**oo | 0 | Because for all complex numbers z near | 110 | | | 0, z**oo -> 0. | 111 +--------------+---------+-----------------------------------------------+ 112 | 0**-oo | zoo | This is not strictly true, as 0**oo may be | 113 | | | oscillating between positive and negative | 114 | | | values or rotating in the complex plane. | 115 | | | It is convenient, however, when the base | 116 | | | is positive. | 117 +--------------+---------+-----------------------------------------------+ 118 | 1**oo | nan | Because there are various cases where | 119 | 1**-oo | | lim(x(t),t)=1, lim(y(t),t)=oo (or -oo), | 120 | | | but lim(x(t)**y(t), t) != 1. See [3]. | 121 +--------------+---------+-----------------------------------------------+ 122 | z**zoo | nan | No limit for z**t for t -> zoo. | 123 +--------------+---------+-----------------------------------------------+ 124 | (-1)**oo | nan | Because of oscillations in the limit. | 125 | (-1)**(-oo) | | | 126 +--------------+---------+-----------------------------------------------+ 127 | oo**oo | oo | | 128 +--------------+---------+-----------------------------------------------+ 129 | oo**-oo | 0 | | 130 +--------------+---------+-----------------------------------------------+ 131 | (-oo)**oo | nan | | 132 | (-oo)**-oo | | | 133 +--------------+---------+-----------------------------------------------+ 134 | oo**I | nan | oo**e could probably be best thought of as | 135 | (-oo)**I | | the limit of x**e for real x as x tends to | 136 | | | oo. If e is I, then the limit does not exist | 137 | | | and nan is used to indicate that. | 138 +--------------+---------+-----------------------------------------------+ 139 | oo**(1+I) | zoo | If the real part of e is positive, then the | 140 | (-oo)**(1+I) | | limit of abs(x**e) is oo. So the limit value | 141 | | | is zoo. | 142 +--------------+---------+-----------------------------------------------+ 143 | oo**(-1+I) | 0 | If the real part of e is negative, then the | 144 | -oo**(-1+I) | | limit is 0. | 145 +--------------+---------+-----------------------------------------------+ 146 147 Because symbolic computations are more flexible that floating point 148 calculations and we prefer to never return an incorrect answer, 149 we choose not to conform to all IEEE 754 conventions. This helps 150 us avoid extra test-case code in the calculation of limits. 151 152 See Also 153 ======== 154 155 diofant.core.numbers.Infinity 156 diofant.core.numbers.NaN 157 158 References 159 ========== 160 161 * https://en.wikipedia.org/wiki/Exponentiation 162 * https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero 163 * https://en.wikipedia.org/wiki/Indeterminate_forms 164 165 """ 166 167 is_Pow = True 168 169 @cacheit 170 def __new__(cls, b, e, evaluate=None): 171 if evaluate is None: 172 evaluate = global_evaluate[0] 173 from ..functions.elementary.exponential import exp_polar 174 175 b = sympify(b, strict=True) 176 e = sympify(e, strict=True) 177 if evaluate: 178 if nan in (b, e): 179 return nan 180 elif e is S.Zero: 181 return Integer(1) 182 elif e is S.One: 183 return b 184 elif e is zoo: 185 return nan 186 elif e.is_integer and _coeff_isneg(b): 187 if e.is_even: 188 b = -b 189 elif e.is_odd: 190 return -Pow(-b, e) 191 if b is S.One: 192 if abs(e).is_infinite: 193 return nan 194 return Integer(1) 195 else: 196 # recognize base as E 197 if not e.is_Atom and b is not E and not isinstance(b, exp_polar): 198 from ..functions import im, log, sign 199 from ..simplify import denom, numer 200 from .exprtools import factor_terms 201 c, ex = factor_terms(e, sign=False).as_coeff_Mul() 202 den = denom(ex) 203 if isinstance(den, log) and den.args[0] == b: 204 return E**(c*numer(ex)) 205 elif den.is_Add: 206 s = sign(im(b)) 207 if s.is_Number and s and den == \ 208 log(-factor_terms(b, sign=False)) + s*I*pi: 209 return E**(c*numer(ex)) 210 211 obj = b._eval_power(e) 212 if obj is not None: 213 return obj 214 return Expr.__new__(cls, b, e) 215 216 def _eval_is_commutative(self): 217 return self.base.is_commutative and self.exp.is_commutative 218 219 @property 220 def base(self): 221 """Returns base of the power expression.""" 222 return self.args[0] 223 224 @property 225 def exp(self): 226 """Returns exponent of the power expression.""" 227 return self.args[1] 228 229 @classmethod 230 def class_key(cls): 231 """Nice order of classes.""" 232 return 4, 2, cls.__name__ 233 234 def _eval_power(self, other): 235 from ..functions import Abs, arg, exp, floor, im, log, re, sign 236 b, e = self.as_base_exp() 237 if b is nan: 238 return (b**e)**other # let __new__ handle it 239 240 s = None 241 if other.is_integer: 242 s = 1 243 elif b.is_polar: # e.g. exp_polar, besselj, var('p', polar=True)... 244 s = 1 245 elif e.is_extended_real is not None: 246 # helper functions =========================== 247 def _half(e): 248 """Return True if the exponent has a literal 2 as the 249 denominator, else None. 250 251 """ 252 if getattr(e, 'denominator', None) == 2: 253 return True 254 n, d = e.as_numer_denom() 255 if n.is_integer and d == 2: 256 return True 257 258 def _n2(e): 259 """Return ``e`` evaluated to a Number with 2 significant 260 digits, else None. 261 262 """ 263 try: 264 rv = e.evalf(2) 265 if rv.is_Number: 266 return rv 267 except PrecisionExhausted: # pragma: no cover 268 pass 269 270 # =================================================== 271 if e.is_extended_real: 272 # we need _half(other) with constant floor or 273 # floor(Rational(1, 2) - e*arg(b)/2/pi) == 0 274 275 # handle -1 as special case 276 if (e == -1): 277 # floor arg. is 1/2 + arg(b)/2/pi 278 if _half(other): 279 if b.is_negative is True: 280 return (-1)**other*Pow(-b, e*other) 281 if b.is_extended_real is False: 282 return Pow(b.conjugate()/Abs(b)**2, other) 283 elif e.is_even: 284 if b.is_extended_real: 285 b = abs(b) 286 if b.is_imaginary: 287 b = abs(im(b))*I 288 289 if (abs(e) < 1) == true or (e == 1): 290 s = 1 # floor = 0 291 elif b.is_nonnegative: 292 s = 1 # floor = 0 293 elif re(b).is_nonnegative and (abs(e) < 2) == true: 294 s = 1 # floor = 0 295 elif im(b).is_nonzero and (abs(e) == 2): 296 s = 1 # floor = 0 297 elif b.is_imaginary and (abs(e) == 2): 298 s = 1 # floor = 0 299 elif _half(other): 300 s = exp(2*pi*I*other*floor( 301 Rational(1, 2) - e*arg(b)/(2*pi))) 302 if s.is_extended_real and _n2(sign(s) - s) == 0: 303 s = sign(s) 304 else: 305 s = None 306 else: 307 # e.is_extended_real is False requires: 308 # _half(other) with constant floor or 309 # floor(Rational(1, 2) - im(e*log(b))/2/pi) == 0 310 s = exp(2*I*pi*other*floor(Rational(1, 2) - im(e*log(b))/2/pi)) 311 # be careful to test that s is -1 or 1 b/c sign(I) == I: 312 # so check that s is real 313 if s.is_extended_real and _n2(sign(s) - s) == 0: 314 s = sign(s) 315 else: 316 s = None 317 318 if s is not None: 319 return s*Pow(b, e*other) 320 321 def _eval_is_positive(self): 322 b, e = self.base, self.exp 323 324 if b.is_nonnegative and b == e: 325 return True 326 elif b.is_positive and (e.is_real or e.is_positive): 327 return True 328 elif b.is_negative and e.is_integer and (b.is_finite or e.is_nonnegative): 329 return e.is_even 330 elif b.is_nonpositive and e.is_odd and (b.is_finite or e.is_nonnegative): 331 return False 332 elif b in {I, -I} and e.is_imaginary: 333 return True 334 335 def _eval_is_nonnegative(self): 336 b, e = self.base, self.exp 337 338 if b.is_imaginary and e.is_nonnegative: 339 m = e % 4 340 if m.is_integer: 341 return m.is_zero 342 343 def _eval_is_negative(self): 344 b, e = self.base, self.exp 345 346 if b.is_negative: 347 if e.is_odd and (b.is_finite or e.is_positive): 348 return True 349 if e.is_even: 350 return False 351 elif b.is_positive: 352 if e.is_extended_real: 353 return False 354 elif b.is_nonnegative: 355 if e.is_nonnegative: 356 return False 357 elif b.is_nonpositive: 358 if e.is_even: 359 return False 360 elif b.is_extended_real: 361 if e.is_even: 362 return False 363 364 def _eval_is_zero(self): 365 b, e = self.base, self.exp 366 367 if b.is_zero: 368 if e.is_positive: 369 return True 370 elif e.is_nonpositive: 371 return False 372 elif b.is_nonzero: 373 if e.is_finite: 374 return False 375 elif e.is_infinite: 376 if (1 - abs(b)).is_positive: 377 return e.is_positive 378 elif (1 - abs(b)).is_negative: 379 return e.is_negative 380 381 def _eval_is_integer(self): 382 b, e = self.base, self.exp 383 384 if b.is_rational: 385 if b.is_integer is False and e.is_positive: 386 return False # rat**nonneg 387 if b.is_integer and e.is_integer: 388 if b is S.NegativeOne: 389 return True 390 if e.is_nonnegative or e.is_positive: 391 return True 392 if b.is_integer and e.is_negative and (e.is_finite or e.is_integer): 393 if (b - 1).is_nonzero and (b + 1).is_nonzero: 394 return False 395 if b.is_Number and e.is_Number: 396 check = self.func(*self.args) 397 if check.is_Integer: 398 return True 399 400 def _eval_is_extended_real(self): 401 from ..functions import arg, log 402 from .mul import Mul 403 404 b, e = self.base, self.exp 405 406 if b is E: 407 if e.is_extended_real: 408 return True 409 elif e.is_imaginary: 410 return (2*I*e/pi).is_even 411 412 if b.is_extended_real is None: 413 if b.func == self.func and b.base is E and b.exp.is_imaginary: 414 return e.is_imaginary 415 if e.is_extended_real is None: 416 return 417 418 if b.is_extended_real and e.is_extended_real: 419 if b.is_positive: 420 return True 421 elif b.is_nonnegative: 422 if e.is_nonnegative: 423 return True 424 else: 425 if e.is_integer: 426 if b.is_nonzero or e.is_nonnegative: 427 return True 428 elif b.is_negative: 429 if e.is_rational and e.is_noninteger: 430 return False 431 432 if b.is_nonzero and e.is_negative: 433 return (b**-e).is_extended_real 434 435 if b.is_imaginary: 436 if e.is_integer: 437 if e.is_even: 438 if b.is_nonzero or e.is_nonnegative: 439 return True 440 elif e.is_odd: 441 return False 442 elif e.is_imaginary and log(b).is_imaginary: 443 return True 444 elif e.is_Add: 445 c, a = e.as_coeff_Add() 446 if c and c.is_Integer: 447 return Mul(b**c, b**a, evaluate=False).is_extended_real 448 elif b in (-I, I) and (e/2).is_noninteger: 449 return False 450 return 451 452 if b.is_extended_real and e.is_imaginary: 453 if b is S.NegativeOne: 454 return True 455 c = e.coeff(I) 456 if c in (1, -1): 457 if b == 2: 458 return False 459 460 if b.is_extended_real is False: # we already know it's not imag 461 i = arg(b)*e/pi 462 return i.is_integer 463 464 def _eval_is_complex(self): 465 from ..functions import log 466 467 b, e = self.base, self.exp 468 469 if b.is_complex: 470 exp = log(b)*e 471 return fuzzy_or([exp.is_complex, exp.is_negative]) 472 473 def _eval_is_imaginary(self): 474 from ..functions import arg, log 475 476 b, e = self.base, self.exp 477 478 if b.is_imaginary: 479 if e.is_integer: 480 return e.is_odd 481 482 if e.is_imaginary and e.is_nonzero: 483 if log(b).is_imaginary: 484 return False 485 486 if b.is_real and e.is_real: 487 if b.is_positive: 488 return False 489 else: 490 if e.is_integer: 491 return False 492 else: 493 if (2*e).is_integer: 494 return b.is_negative 495 496 if b.is_real is False: # we already know it's not imag 497 return (2*arg(b)*e/pi).is_odd 498 499 def _eval_is_odd(self): 500 b, e = self.base, self.exp 501 502 if e.is_integer: 503 if e.is_positive: 504 return b.is_odd 505 elif e.is_nonnegative and b.is_odd: 506 return True 507 elif b is S.NegativeOne: 508 return True 509 510 def _eval_is_finite(self): 511 b, e = self.base, self.exp 512 513 if e.is_negative: 514 if b.is_zero: 515 return False 516 if b.is_finite and e.is_finite: 517 if e.is_nonnegative or b.is_nonzero: 518 return True 519 520 def _eval_is_polar(self): 521 b, e = self.base, self.exp 522 523 if b.is_polar and e.is_commutative: 524 return True 525 526 def _eval_subs(self, old, new): 527 from ..functions import log 528 from .symbol import Symbol 529 530 def _check(ct1, ct2, old): 531 """Return bool, pow where, if bool is True, then the exponent of 532 Pow `old` will combine with `pow` so the substitution is valid, 533 otherwise bool will be False, 534 535 cti are the coefficient and terms of an exponent of self or old 536 In this _eval_subs routine a change like (b**(2*x)).subs({b**x: y}) 537 will give y**2 since (b**x)**2 == b**(2*x); if that equality does 538 not hold then the substitution should not occur so `bool` will be 539 False. 540 541 """ 542 coeff1, terms1 = ct1 543 coeff2, terms2 = ct2 544 if terms1 == terms2: 545 pow = coeff1/coeff2 546 try: 547 pow = as_int(pow) 548 combines = True 549 except ValueError: 550 combines = Pow._eval_power( 551 Pow(*old.as_base_exp(), evaluate=False), 552 pow) 553 if isinstance(combines, Pow): 554 combines = combines.base is old.base 555 else: 556 combines = False 557 return combines, pow 558 return False, None 559 560 if old == self.base: 561 return new**self.exp._subs(old, new) 562 563 if old.func is self.func and self.exp == old.exp: 564 l = log(self.base, old.base) 565 if l.is_Number: 566 return Pow(new, l) 567 568 if isinstance(old, self.func) and self.base == old.base: 569 if self.exp.is_Add is False: 570 ct2 = old.exp.as_independent(Symbol, as_Add=False) 571 ct1 = (self.exp/ct2[1], ct2[1]) 572 ok, pow = _check(ct1, ct2, old) 573 if ok: 574 # issue sympy/sympy#5180: (x**(6*y)).subs({x**(3*y):z})->z**2 575 return self.func(new, pow) 576 else: # b**(6*x+a).subs({b**(3*x): y}) -> y**2 * b**a 577 # exp(exp(x) + exp(x**2)).subs({exp(exp(x)): w}) -> w * exp(exp(x**2)) 578 oarg = old.exp 579 new_l = [] 580 o_al = [] 581 ct2 = oarg.as_coeff_mul() 582 for a in self.exp.args: 583 newa = a._subs(old, new) 584 ct1 = newa.as_coeff_mul() 585 ok, pow = _check(ct1, ct2, old) 586 if ok: 587 new_l.append(new**pow) 588 continue 589 o_al.append(newa) 590 if new_l: 591 new_l.append(Pow(self.base, Add(*o_al), evaluate=False)) 592 return Mul(*new_l) 593 594 if old.is_Pow and old.base is E and self.exp.is_extended_real and self.base.is_positive: 595 ct1 = old.exp.as_independent(Symbol, as_Add=False) 596 ct2 = (self.exp*log(self.base)).as_independent( 597 Symbol, as_Add=False) 598 ok, pow = _check(ct1, ct2, old) 599 if ok: 600 return self.func(new, pow) # (2**x).subs({exp(x*log(2)): z}) -> z 601 602 def as_base_exp(self): 603 """Return base and exp of self. 604 605 If base is 1/Integer, then return Integer, -exp. If this extra 606 processing is not needed, the base and exp properties will 607 give the raw arguments 608 609 Examples 610 ======== 611 612 >>> p = Pow(Rational(1, 2), 2, evaluate=False) 613 >>> p.as_base_exp() 614 (2, -2) 615 >>> p.args 616 (1/2, 2) 617 618 """ 619 b, e = self.base, self.exp 620 if b.is_Rational and b.numerator == 1 and b.denominator != 1: 621 return Integer(b.denominator), -e 622 return b, e 623 624 def _eval_adjoint(self): 625 from ..functions.elementary.complexes import adjoint 626 i, p = self.exp.is_integer, self.base.is_positive 627 if i: 628 return adjoint(self.base)**self.exp 629 if p: 630 return self.base**adjoint(self.exp) 631 632 def _eval_conjugate(self): 633 if self.is_extended_real: 634 return self 635 from ..functions.elementary.complexes import conjugate as c 636 i, p = self.exp.is_integer, self.base.is_positive 637 if i: 638 return c(self.base)**self.exp 639 if p: 640 return self.base**c(self.exp) 641 if i is False and p is False: 642 expanded = expand_complex(self) 643 assert expanded != self 644 return c(expanded) 645 646 def _eval_transpose(self): 647 from ..functions.elementary.complexes import transpose 648 i, p = self.exp.is_integer, self.base.is_complex 649 if p: 650 return self.base**self.exp 651 if i: 652 return transpose(self.base)**self.exp 653 654 def _eval_expand_power_exp(self, **hints): 655 """a**(n+m) -> a**n*a**m.""" 656 b = self.base 657 e = self.exp 658 if e.is_Add and e.is_commutative: 659 expr = [] 660 for x in e.args: 661 expr.append(self.func(self.base, x)) 662 return Mul(*expr) 663 return self.func(b, e) 664 665 def _eval_expand_power_base(self, **hints): 666 """(a*b)**n -> a**n * b**n.""" 667 force = hints.get('force', False) 668 669 b = self.base 670 e = self.exp 671 if not b.is_Mul: 672 return self 673 674 cargs, nc = b.args_cnc(split_1=False) 675 676 # expand each term - this is top-level-only 677 # expansion but we have to watch out for things 678 # that don't have an _eval_expand method 679 if nc: 680 nc = [i._eval_expand_power_base(**hints) 681 if hasattr(i, '_eval_expand_power_base') else i 682 for i in nc] 683 684 if e.is_Integer: 685 if e.is_positive: 686 rv = Mul(*nc*e) 687 else: 688 rv = 1/Mul(*nc*-e) 689 if cargs: 690 rv *= Mul(*cargs)**e 691 return rv 692 693 if not cargs: 694 return self.func(Mul(*nc), e, evaluate=False) 695 696 nc = [Mul(*nc)] 697 698 # sift the commutative bases 699 def pred(x): 700 if x is I: 701 return I 702 polar = x.is_polar 703 if polar: 704 return True 705 if polar is None: 706 return fuzzy_or([x.is_nonnegative, (1/x).is_nonnegative]) 707 sifted = sift(cargs, pred) 708 nonneg = sifted[True] 709 other = sifted[None] 710 neg = sifted[False] 711 imag = sifted[I] 712 if imag: 713 i = len(imag) % 4 714 if i == 0: 715 pass 716 elif i == 1: 717 other.append(I) 718 elif i == 2: 719 if neg: 720 nonn = -neg.pop() 721 if nonn is not S.One: 722 nonneg.append(nonn) 723 else: 724 neg.append(Integer(-1)) 725 else: 726 if neg: 727 nonn = -neg.pop() 728 if nonn is not S.One: 729 nonneg.append(nonn) 730 else: 731 neg.append(Integer(-1)) 732 other.append(I) 733 del imag 734 735 # bring out the bases that can be separated from the base 736 737 if force or e.is_integer: 738 # treat all commutatives the same and put nc in other 739 cargs = nonneg + neg + other 740 other = nc 741 else: 742 # this is just like what is happening automatically, except 743 # that now we are doing it for an arbitrary exponent for which 744 # no automatic expansion is done 745 746 assert not e.is_Integer 747 748 # handle negatives by making them all positive and putting 749 # the residual -1 in other 750 if len(neg) > 1: 751 o = Integer(1) 752 if not other and neg[0].is_Number: 753 o *= neg.pop(0) 754 if len(neg) % 2: 755 o = -o 756 for n in neg: 757 nonneg.append(-n) 758 if o is not S.One: 759 other.append(o) 760 elif neg and other: 761 if neg[0].is_Number and neg[0] is not S.NegativeOne: 762 other.append(Integer(-1)) 763 nonneg.append(-neg[0]) 764 else: 765 other.extend(neg) 766 else: 767 other.extend(neg) 768 del neg 769 770 cargs = nonneg 771 other += nc 772 773 rv = Integer(1) 774 if cargs: 775 rv *= Mul(*[self.func(b, e, evaluate=False) for b in cargs]) 776 if other: 777 rv *= self.func(Mul(*other), e, evaluate=False) 778 return rv 779 780 def _eval_expand_multinomial(self, **hints): 781 """(a+b+..) ** n -> a**n + n*a**(n-1)*b + .., n is nonzero integer.""" 782 base, exp = self.base, self.exp 783 result = self 784 785 if exp.is_Rational and exp.numerator > 0 and base.is_Add: 786 if not exp.is_Integer: 787 n = Integer(exp.numerator // exp.denominator) 788 789 if not n: 790 return result 791 else: 792 radical, result = self.func(base, exp - n), [] 793 794 expanded_base_n = self.func(base, n) 795 if expanded_base_n.is_Pow: 796 expanded_base_n = \ 797 expanded_base_n._eval_expand_multinomial() 798 for term in Add.make_args(expanded_base_n): 799 result.append(term*radical) 800 801 return Add(*result) 802 803 n = int(exp) 804 805 if base.is_commutative: 806 order_terms, other_terms = [], [] 807 808 for b in base.args: 809 if b.is_Order: 810 order_terms.append(b) 811 else: 812 other_terms.append(b) 813 814 if order_terms: 815 # (f(x) + O(x^n))^m -> f(x)^m + m*f(x)^{m-1} *O(x^n) 816 f = Add(*other_terms) 817 o = Add(*order_terms) 818 819 if n == 2: 820 return expand_multinomial(f**n, deep=False) + n*f*o 821 else: 822 g = expand_multinomial(f**(n - 1), deep=False) 823 return expand_mul(f*g, deep=False) + n*g*o 824 825 if base.is_number: 826 # Efficiently expand expressions of the form (a + b*I)**n 827 # where 'a' and 'b' are real numbers and 'n' is integer. 828 a, b = base.as_real_imag() 829 830 if a.is_Rational and b.is_Rational: 831 if not a.is_Integer: 832 if not b.is_Integer: 833 k = self.func(a.denominator * b.denominator, n) 834 a, b = a.numerator*b.denominator, a.denominator*b.numerator 835 else: 836 k = self.func(a.denominator, n) 837 a, b = a.numerator, a.denominator*b 838 elif not b.is_Integer: 839 k = self.func(b.denominator, n) 840 a, b = a*b.denominator, b.numerator 841 else: 842 k = 1 843 844 a, b, c, d = int(a), int(b), 1, 0 845 846 while n: 847 if n & 1: 848 c, d = a*c - b*d, b*c + a*d 849 n -= 1 850 a, b = a*a - b*b, 2*a*b 851 n //= 2 852 853 if k == 1: 854 return c + I*d 855 else: 856 return Integer(c)/k + I*d/k 857 858 p = other_terms 859 # (x+y)**3 -> x**3 + 3*x**2*y + 3*x*y**2 + y**3 860 # in this particular example: 861 # p = [x,y]; n = 3 862 # so now it's easy to get the correct result -- we get the 863 # coefficients first: 864 from ..ntheory import multinomial_coefficients 865 from ..polys import Poly 866 expansion_dict = multinomial_coefficients(len(p), n) 867 # in our example: {(3, 0): 1, (1, 2): 3, (0, 3): 1, (2, 1): 3} 868 # and now construct the expression. 869 return Poly(expansion_dict, *p).as_expr() 870 else: 871 if n == 2: 872 return Add(*[f*g for f in base.args for g in base.args]) 873 else: 874 multi = (base**(n - 1))._eval_expand_multinomial() 875 assert multi.is_Add 876 return Add(*[f*g for f in base.args for g in multi.args]) 877 elif (exp.is_Rational and exp.numerator < 0 and base.is_Add and 878 abs(exp.numerator) > exp.denominator): 879 return 1 / self.func(base, -exp)._eval_expand_multinomial() 880 elif exp.is_Add and base.is_Number: 881 # a + b a b 882 # n --> n n , where n, a, b are Numbers 883 884 coeff, tail = Integer(1), Integer(0) 885 for term in exp.args: 886 if term.is_Number: 887 coeff *= self.func(base, term) 888 else: 889 tail += term 890 891 return coeff * self.func(base, tail) 892 else: 893 return result 894 895 def as_real_imag(self, deep=True, **hints): 896 """Returns real and imaginary parts of self 897 898 See Also 899 ======== 900 901 diofant.core.expr.Expr.as_real_imag 902 903 """ 904 from ..functions import arg, cos, sin 905 906 if self.exp.is_Integer: 907 exp = self.exp 908 re, im = self.base.as_real_imag(deep=deep) 909 if not im: 910 return self, Integer(0) 911 a, b = symbols('a b', cls=Dummy) 912 if exp >= 0: 913 if re.is_Number and im.is_Number: 914 # We can be more efficient in this case 915 expr = expand_multinomial(self.base**exp) 916 return expr.as_real_imag() 917 918 expr = ((a + b)**exp).as_poly() # a = re, b = im; expr = (a + b*I)**exp 919 else: 920 mag = re**2 + im**2 921 re, im = re/mag, -im/mag 922 if re.is_Number and im.is_Number: 923 # We can be more efficient in this case 924 expr = expand_multinomial((re + im*I)**-exp) 925 return expr.as_real_imag() 926 927 expr = ((a + b)**-exp).as_poly() 928 929 # Terms with even b powers will be real 930 r = [i for i in expr.terms() if not i[0][1] % 2] 931 re_part = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r]) 932 # Terms with odd b powers will be imaginary 933 r = [i for i in expr.terms() if i[0][1] % 4 == 1] 934 im_part1 = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r]) 935 r = [i for i in expr.terms() if i[0][1] % 4 == 3] 936 im_part3 = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r]) 937 938 return (re_part.subs({a: re, b: I*im}), 939 im_part1.subs({a: re, b: im}) + im_part3.subs({a: re, b: -im})) 940 941 elif self.exp.is_Rational: 942 re, im = self.base.as_real_imag(deep=deep) 943 944 if im.is_zero and self.exp is S.Half: 945 if re.is_nonnegative: 946 return self, Integer(0) 947 if re.is_nonpositive: 948 return Integer(0), (-self.base)**self.exp 949 950 # XXX: This is not totally correct since for x**(p/q) with 951 # x being imaginary there are actually q roots, but 952 # only a single one is returned from here. 953 r = self.func(self.func(re, 2) + self.func(im, 2), Rational(1, 2)) 954 t = arg(re + I*im) 955 956 rp, tp = self.func(r, self.exp), t*self.exp 957 958 return rp*cos(tp), rp*sin(tp) 959 elif self.base is E: 960 from ..functions import exp 961 re, im = self.exp.as_real_imag() 962 if deep: 963 re = re.expand(deep, **hints) 964 im = im.expand(deep, **hints) 965 c, s = cos(im), sin(im) 966 return exp(re)*c, exp(re)*s 967 else: 968 from ..functions import im, re 969 if deep: 970 hints['complex'] = False 971 972 expanded = self.expand(deep, **hints) 973 if hints.get('ignore') != expanded: 974 return re(expanded), im(expanded) 975 else: 976 return re(self), im(self) 977 978 def _eval_derivative(self, s): 979 from ..functions import log 980 dbase = self.base.diff(s) 981 dexp = self.exp.diff(s) 982 return self * (dexp * log(self.base) + dbase * self.exp/self.base) 983 984 def _eval_evalf(self, prec): 985 base, exp = self.as_base_exp() 986 base = base._evalf(prec) 987 if not exp.is_Integer: 988 exp = exp._evalf(prec) 989 return self.func(base, exp) 990 991 def _eval_is_polynomial(self, syms): 992 if self.exp.has(*syms): 993 return False 994 995 if self.base.has(*syms): 996 return bool(self.base._eval_is_polynomial(syms) and 997 self.exp.is_Integer and (self.exp >= 0)) 998 else: 999 return True 1000 1001 def _eval_is_rational(self): 1002 p = self.func(*self.as_base_exp()) # in case it's unevaluated 1003 if not p.is_Pow: 1004 return p.is_rational 1005 b, e = p.base, p.exp 1006 1007 if e.is_Rational and b.is_Rational: 1008 # we didn't check that e is not an Integer 1009 # because Rational**Integer autosimplifies 1010 return False 1011 if e.is_integer: 1012 if b.is_rational: 1013 if b.is_nonzero or e.is_nonnegative: 1014 return True 1015 if b == e: # always rational, even for 0**0 1016 return True 1017 elif b.is_irrational: 1018 if e.is_zero: 1019 return True 1020 if b is E: 1021 if e.is_rational and e.is_nonzero: 1022 return False 1023 1024 def _eval_is_algebraic(self): 1025 b, e = self.base, self.exp 1026 1027 if b.is_zero or (b - 1).is_zero: 1028 return True 1029 elif b is E: 1030 s = self.doit() 1031 if s.func == self.func: 1032 if e.is_nonzero: 1033 if e.is_algebraic: 1034 return False 1035 elif (e/pi).is_rational: 1036 return False 1037 elif (e/(I*pi)).is_rational: 1038 return True 1039 else: 1040 return s.is_algebraic 1041 elif e.is_rational and e.is_nonzero: 1042 if b.is_nonzero or e.is_nonnegative: 1043 return b.is_algebraic 1044 elif b.is_algebraic and e.is_algebraic: 1045 if (b.is_nonzero and (b - 1).is_nonzero) or b.is_irrational: 1046 return e.is_rational 1047 1048 def _eval_is_rational_function(self, syms): 1049 if self.exp.has(*syms): 1050 return False 1051 1052 if self.base.has(*syms): 1053 return self.base._eval_is_rational_function(syms) and \ 1054 self.exp.is_Integer 1055 else: 1056 return True 1057 1058 def _eval_is_algebraic_expr(self, syms): 1059 if self.exp.has(*syms): 1060 return False 1061 1062 if self.base.has(*syms): 1063 return self.base._eval_is_algebraic_expr(syms) and \ 1064 self.exp.is_Rational 1065 else: 1066 return True 1067 1068 def _eval_as_numer_denom(self): 1069 """Expression -> a/b -> a, b. 1070 1071 See Also 1072 ======== 1073 1074 diofant.core.expr.Expr.as_numer_denom 1075 1076 """ 1077 if not self.is_commutative: 1078 return self, Integer(1) 1079 base, exp = self.as_base_exp() 1080 if base is S.One: 1081 return self, Integer(1) 1082 n, d = base.as_numer_denom() 1083 neg_exp = exp.is_negative 1084 if not neg_exp and not (-exp).is_negative: 1085 neg_exp = _coeff_isneg(exp) 1086 int_exp = exp.is_integer 1087 # the denominator cannot be separated from the numerator if 1088 # its sign is unknown unless the exponent is an integer, e.g. 1089 # sqrt(a/b) != sqrt(a)/sqrt(b) when a=1 and b=-1. But if the 1090 # denominator is negative the numerator and denominator can 1091 # be negated and the denominator (now positive) separated. 1092 if not (d.is_extended_real or int_exp): 1093 n = base 1094 d = Integer(1) 1095 dnonpos = d.is_nonpositive 1096 if dnonpos: 1097 n, d = -n, -d 1098 elif dnonpos is None and not int_exp: 1099 n = base 1100 d = Integer(1) 1101 if neg_exp: 1102 n, d = d, n 1103 exp = -exp 1104 if d is S.One: 1105 return self.func(n, exp), Integer(1) 1106 if n is S.One: 1107 return Integer(1), self.func(d, exp) 1108 return self.func(n, exp), self.func(d, exp) 1109 1110 def _matches(self, expr, repl_dict={}): 1111 """Helper method for match(). 1112 1113 See Also 1114 ======== 1115 1116 diofant.core.basic.Basic.matches 1117 1118 """ 1119 expr = sympify(expr, strict=True) 1120 1121 # special case, pattern = 1 and expr.exp can match to 0 1122 if expr is S.One: 1123 d = repl_dict.copy() 1124 d = self.exp._matches(Integer(0), d) 1125 if d is not None: 1126 return d 1127 1128 # make sure the expression to be matched is an Expr 1129 if not isinstance(expr, Expr): 1130 return 1131 1132 b, e = expr.as_base_exp() 1133 1134 # special case number 1135 sb, se = self.as_base_exp() 1136 if sb.is_Symbol and se.is_Integer and expr: 1137 if e.is_rational: 1138 return sb._matches(b**(e/se), repl_dict) 1139 return sb._matches(expr**(1/se), repl_dict) 1140 1141 d = repl_dict.copy() 1142 d = self.base._matches(b, d) 1143 if d is None: 1144 return 1145 1146 d = self.exp.xreplace(d)._matches(e, d) 1147 if d is None: 1148 return Expr._matches(self, expr, repl_dict) 1149 return d 1150 1151 def _eval_nseries(self, x, n, logx): 1152 from ..functions import arg, exp, floor, log 1153 from ..series import Order, limit 1154 from ..simplify import powsimp 1155 if self.base is E: 1156 e_series = self.exp.nseries(x, n=n, logx=logx) 1157 if e_series.is_Order: 1158 return 1 + e_series 1159 e0 = limit(e_series.removeO(), x, 0) 1160 if e0 in (-oo, oo): 1161 return self 1162 t = e_series - e0 1163 exp_series = term = exp(e0) 1164 # series of exp(e0 + t) in t 1165 for i in range(1, n): 1166 term *= t/i 1167 term = term.nseries(x, n=n, logx=logx) 1168 exp_series += term 1169 exp_series += Order(t**n, x) 1170 return powsimp(exp_series, deep=True, combine='exp') 1171 elif self.exp.has(x): 1172 return exp(self.exp*log(self.base)).nseries(x, n=n, logx=logx) 1173 else: 1174 b_series = self.base.nseries(x, n=n, logx=logx) 1175 while b_series.is_Order: 1176 n += 1 1177 b_series = self.base.nseries(x, n=n, logx=logx) 1178 b0 = b_series.as_leading_term(x) 1179 t = expand_mul((b_series/b0 - 1).cancel()) 1180 if t.is_Add: 1181 t = t.func(*[i for i in t.args if i.limit(x, 0).is_finite]) 1182 c, e = b0.as_coeff_exponent(x) 1183 if self.exp is oo: 1184 if e != 0: 1185 sig = -e 1186 else: 1187 sig = abs(c) - 1 if c != 1 else t.removeO() 1188 if sig.is_positive: 1189 return oo 1190 elif sig.is_negative: 1191 return Integer(0) 1192 else: 1193 raise NotImplementedError 1194 pow_series = term = Integer(1) 1195 # series of (1 + t)**e in t 1196 for i in range(1, n): 1197 term *= (self.exp - i + 1)*t/i 1198 term = term.nseries(x, n=n, logx=logx) 1199 pow_series += term 1200 factor = b0**self.exp 1201 if t != 0 and not (self.exp.is_Integer and self.exp >= 0 and n > self.exp): 1202 pow_series += Order(t**n, x) 1203 # branch handling 1204 if c.is_negative: 1205 l = floor(arg(t.removeO()*c)/(2*pi)).limit(x, 0) 1206 assert l.is_finite 1207 factor *= exp(2*pi*I*self.exp*l) 1208 pow_series = expand_mul(factor*pow_series) 1209 return powsimp(pow_series, deep=True, combine='exp') 1210 1211 def _eval_as_leading_term(self, x): 1212 from ..functions import exp, log 1213 from ..series import Order 1214 if not self.exp.has(x): 1215 return self.func(self.base.as_leading_term(x), self.exp) 1216 elif self.base is E: 1217 if self.exp.is_Mul: 1218 k, arg = self.exp.as_independent(x) 1219 else: 1220 k, arg = Integer(1), self.exp 1221 if arg.is_Add: 1222 return Mul(*[exp(k*f).as_leading_term(x) for f in arg.args]) 1223 arg = self.exp.as_leading_term(x) 1224 if Order(1, x).contains(arg): 1225 return Integer(1) 1226 return exp(arg) 1227 else: 1228 return exp(self.exp*log(self.base)).as_leading_term(x) 1229 1230 def _eval_rewrite_as_sin(self, base, exp): 1231 from ..functions import sin 1232 if self.base is E: 1233 return sin(I*self.exp + pi/2) - I*sin(I*self.exp) 1234 1235 def _eval_rewrite_as_cos(self, base, exp): 1236 from ..functions import cos 1237 if self.base is E: 1238 return cos(I*self.exp) + I*cos(I*self.exp + pi/2) 1239 1240 def _eval_rewrite_as_tanh(self, base, exp): 1241 from ..functions import tanh 1242 if self.base is E: 1243 return (1 + tanh(self.exp/2))/(1 - tanh(self.exp/2)) 1244 1245 def as_content_primitive(self, radical=False): 1246 """Return the tuple (R, self/R) where R is the positive Rational 1247 extracted from self. 1248 1249 Examples 1250 ======== 1251 1252 >>> sqrt(4 + 4*sqrt(2)).as_content_primitive() 1253 (2, sqrt(1 + sqrt(2))) 1254 >>> sqrt(3 + 3*sqrt(2)).as_content_primitive() 1255 (1, sqrt(3)*sqrt(1 + sqrt(2))) 1256 1257 >>> ((2*x + 2)**2).as_content_primitive() 1258 (4, (x + 1)**2) 1259 >>> (4**((1 + y)/2)).as_content_primitive() 1260 (2, 4**(y/2)) 1261 >>> (3**((1 + y)/2)).as_content_primitive() 1262 (1, 3**((y + 1)/2)) 1263 >>> (3**((5 + y)/2)).as_content_primitive() 1264 (9, 3**((y + 1)/2)) 1265 >>> eq = 3**(2 + 2*x) 1266 >>> powsimp(eq) == eq 1267 True 1268 >>> eq.as_content_primitive() 1269 (9, 3**(2*x)) 1270 >>> powsimp(Mul(*_)) 1271 3**(2*x + 2) 1272 1273 >>> eq = (2 + 2*x)**y 1274 >>> s = expand_power_base(eq) 1275 >>> s.is_Mul, s 1276 (False, (2*x + 2)**y) 1277 >>> eq.as_content_primitive() 1278 (1, (2*(x + 1))**y) 1279 >>> s = expand_power_base(_[1]) 1280 >>> s.is_Mul, s 1281 (True, 2**y*(x + 1)**y) 1282 1283 See Also 1284 ======== 1285 1286 diofant.core.expr.Expr.as_content_primitive 1287 1288 """ 1289 b, e = self.as_base_exp() 1290 b = _keep_coeff(*b.as_content_primitive(radical=radical)) 1291 ce, pe = e.as_content_primitive(radical=radical) 1292 if b.is_Rational: 1293 # e 1294 # = ce*pe 1295 # = ce*(h + t) 1296 # = ce*h + ce*t 1297 # => self 1298 # = b**(ce*h)*b**(ce*t) 1299 # = b**(cehp/cehq)*b**(ce*t) 1300 # = b**(iceh+r/cehq)*b**(ce*t) 1301 # = b**iceh*b**(r/cehq)*b**(ce*t) 1302 # = b**iceh*b**(ce*t + r/cehq) 1303 h, t = pe.as_coeff_Add() 1304 if h.is_Rational: 1305 ceh = ce*h 1306 c = self.func(b, ceh) 1307 r = Integer(0) 1308 if not c.is_Rational: 1309 iceh, r = divmod(ceh.numerator, ceh.denominator) 1310 c = self.func(b, iceh) 1311 return c, self.func(b, _keep_coeff(ce, t + r/ce/ceh.denominator)) 1312 e = _keep_coeff(ce, pe) 1313 # b**e = (h*t)**e = h**e*t**e = c*m*t**e 1314 if e.is_Rational and b.is_Mul: 1315 h, t = b.as_content_primitive(radical=radical) # h is positive 1316 c, m = self.func(h, e).as_coeff_Mul() # so c is positive 1317 m, me = m.as_base_exp() 1318 if m is S.One or me == e: # probably always true 1319 # return the following, not return c, m*Pow(t, e) 1320 # which would change Pow into Mul; we let diofant 1321 # decide what to do by using the unevaluated Mul, e.g 1322 # should it stay as sqrt(2 + 2*sqrt(5)) or become 1323 # sqrt(2)*sqrt(1 + sqrt(5)) 1324 return c, self.func(_keep_coeff(m, t), e) 1325 return Integer(1), self.func(b, e) 1326