Coverage for /usr/lib/python3/dist-packages/sympy/core/power.py: 26%
1321 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
2from typing import Callable
3from math import log as _log, sqrt as _sqrt
4from itertools import product
6from .sympify import _sympify
7from .cache import cacheit
8from .singleton import S
9from .expr import Expr
10from .evalf import PrecisionExhausted
11from .function import (expand_complex, expand_multinomial,
12 expand_mul, _mexpand, PoleError)
13from .logic import fuzzy_bool, fuzzy_not, fuzzy_and, fuzzy_or
14from .parameters import global_parameters
15from .relational import is_gt, is_lt
16from .kind import NumberKind, UndefinedKind
17from sympy.external.gmpy import HAS_GMPY, gmpy
18from sympy.utilities.iterables import sift
19from sympy.utilities.exceptions import sympy_deprecation_warning
20from sympy.utilities.misc import as_int
21from sympy.multipledispatch import Dispatcher
23from mpmath.libmp import sqrtrem as mpmath_sqrtrem
27def isqrt(n):
28 """Return the largest integer less than or equal to sqrt(n)."""
29 if n < 0:
30 raise ValueError("n must be nonnegative")
31 n = int(n)
33 # Fast path: with IEEE 754 binary64 floats and a correctly-rounded
34 # math.sqrt, int(math.sqrt(n)) works for any integer n satisfying 0 <= n <
35 # 4503599761588224 = 2**52 + 2**27. But Python doesn't guarantee either
36 # IEEE 754 format floats *or* correct rounding of math.sqrt, so check the
37 # answer and fall back to the slow method if necessary.
38 if n < 4503599761588224:
39 s = int(_sqrt(n))
40 if 0 <= n - s*s <= 2*s:
41 return s
43 return integer_nthroot(n, 2)[0]
46def integer_nthroot(y, n):
47 """
48 Return a tuple containing x = floor(y**(1/n))
49 and a boolean indicating whether the result is exact (that is,
50 whether x**n == y).
52 Examples
53 ========
55 >>> from sympy import integer_nthroot
56 >>> integer_nthroot(16, 2)
57 (4, True)
58 >>> integer_nthroot(26, 2)
59 (5, False)
61 To simply determine if a number is a perfect square, the is_square
62 function should be used:
64 >>> from sympy.ntheory.primetest import is_square
65 >>> is_square(26)
66 False
68 See Also
69 ========
70 sympy.ntheory.primetest.is_square
71 integer_log
72 """
73 y, n = as_int(y), as_int(n)
74 if y < 0:
75 raise ValueError("y must be nonnegative")
76 if n < 1:
77 raise ValueError("n must be positive")
78 if HAS_GMPY and n < 2**63:
79 # Currently it works only for n < 2**63, else it produces TypeError
80 # sympy issue: https://github.com/sympy/sympy/issues/18374
81 # gmpy2 issue: https://github.com/aleaxit/gmpy/issues/257
82 if HAS_GMPY >= 2:
83 x, t = gmpy.iroot(y, n)
84 else:
85 x, t = gmpy.root(y, n)
86 return as_int(x), bool(t)
87 return _integer_nthroot_python(y, n)
89def _integer_nthroot_python(y, n):
90 if y in (0, 1):
91 return y, True
92 if n == 1:
93 return y, True
94 if n == 2:
95 x, rem = mpmath_sqrtrem(y)
96 return int(x), not rem
97 if n >= y.bit_length():
98 return 1, False
99 # Get initial estimate for Newton's method. Care must be taken to
100 # avoid overflow
101 try:
102 guess = int(y**(1./n) + 0.5)
103 except OverflowError:
104 exp = _log(y, 2)/n
105 if exp > 53:
106 shift = int(exp - 53)
107 guess = int(2.0**(exp - shift) + 1) << shift
108 else:
109 guess = int(2.0**exp)
110 if guess > 2**50:
111 # Newton iteration
112 xprev, x = -1, guess
113 while 1:
114 t = x**(n - 1)
115 xprev, x = x, ((n - 1)*x + y//t)//n
116 if abs(x - xprev) < 2:
117 break
118 else:
119 x = guess
120 # Compensate
121 t = x**n
122 while t < y:
123 x += 1
124 t = x**n
125 while t > y:
126 x -= 1
127 t = x**n
128 return int(x), t == y # int converts long to int if possible
131def integer_log(y, x):
132 r"""
133 Returns ``(e, bool)`` where e is the largest nonnegative integer
134 such that :math:`|y| \geq |x^e|` and ``bool`` is True if $y = x^e$.
136 Examples
137 ========
139 >>> from sympy import integer_log
140 >>> integer_log(125, 5)
141 (3, True)
142 >>> integer_log(17, 9)
143 (1, False)
144 >>> integer_log(4, -2)
145 (2, True)
146 >>> integer_log(-125,-5)
147 (3, True)
149 See Also
150 ========
151 integer_nthroot
152 sympy.ntheory.primetest.is_square
153 sympy.ntheory.factor_.multiplicity
154 sympy.ntheory.factor_.perfect_power
155 """
156 if x == 1:
157 raise ValueError('x cannot take value as 1')
158 if y == 0:
159 raise ValueError('y cannot take value as 0')
161 if x in (-2, 2):
162 x = int(x)
163 y = as_int(y)
164 e = y.bit_length() - 1
165 return e, x**e == y
166 if x < 0:
167 n, b = integer_log(y if y > 0 else -y, -x)
168 return n, b and bool(n % 2 if y < 0 else not n % 2)
170 x = as_int(x)
171 y = as_int(y)
172 r = e = 0
173 while y >= x:
174 d = x
175 m = 1
176 while y >= d:
177 y, rem = divmod(y, d)
178 r = r or rem
179 e += m
180 if y > d:
181 d *= d
182 m *= 2
183 return e, r == 0 and y == 1
186class Pow(Expr):
187 """
188 Defines the expression x**y as "x raised to a power y"
190 .. deprecated:: 1.7
192 Using arguments that aren't subclasses of :class:`~.Expr` in core
193 operators (:class:`~.Mul`, :class:`~.Add`, and :class:`~.Pow`) is
194 deprecated. See :ref:`non-expr-args-deprecated` for details.
196 Singleton definitions involving (0, 1, -1, oo, -oo, I, -I):
198 +--------------+---------+-----------------------------------------------+
199 | expr | value | reason |
200 +==============+=========+===============================================+
201 | z**0 | 1 | Although arguments over 0**0 exist, see [2]. |
202 +--------------+---------+-----------------------------------------------+
203 | z**1 | z | |
204 +--------------+---------+-----------------------------------------------+
205 | (-oo)**(-1) | 0 | |
206 +--------------+---------+-----------------------------------------------+
207 | (-1)**-1 | -1 | |
208 +--------------+---------+-----------------------------------------------+
209 | S.Zero**-1 | zoo | This is not strictly true, as 0**-1 may be |
210 | | | undefined, but is convenient in some contexts |
211 | | | where the base is assumed to be positive. |
212 +--------------+---------+-----------------------------------------------+
213 | 1**-1 | 1 | |
214 +--------------+---------+-----------------------------------------------+
215 | oo**-1 | 0 | |
216 +--------------+---------+-----------------------------------------------+
217 | 0**oo | 0 | Because for all complex numbers z near |
218 | | | 0, z**oo -> 0. |
219 +--------------+---------+-----------------------------------------------+
220 | 0**-oo | zoo | This is not strictly true, as 0**oo may be |
221 | | | oscillating between positive and negative |
222 | | | values or rotating in the complex plane. |
223 | | | It is convenient, however, when the base |
224 | | | is positive. |
225 +--------------+---------+-----------------------------------------------+
226 | 1**oo | nan | Because there are various cases where |
227 | 1**-oo | | lim(x(t),t)=1, lim(y(t),t)=oo (or -oo), |
228 | | | but lim( x(t)**y(t), t) != 1. See [3]. |
229 +--------------+---------+-----------------------------------------------+
230 | b**zoo | nan | Because b**z has no limit as z -> zoo |
231 +--------------+---------+-----------------------------------------------+
232 | (-1)**oo | nan | Because of oscillations in the limit. |
233 | (-1)**(-oo) | | |
234 +--------------+---------+-----------------------------------------------+
235 | oo**oo | oo | |
236 +--------------+---------+-----------------------------------------------+
237 | oo**-oo | 0 | |
238 +--------------+---------+-----------------------------------------------+
239 | (-oo)**oo | nan | |
240 | (-oo)**-oo | | |
241 +--------------+---------+-----------------------------------------------+
242 | oo**I | nan | oo**e could probably be best thought of as |
243 | (-oo)**I | | the limit of x**e for real x as x tends to |
244 | | | oo. If e is I, then the limit does not exist |
245 | | | and nan is used to indicate that. |
246 +--------------+---------+-----------------------------------------------+
247 | oo**(1+I) | zoo | If the real part of e is positive, then the |
248 | (-oo)**(1+I) | | limit of abs(x**e) is oo. So the limit value |
249 | | | is zoo. |
250 +--------------+---------+-----------------------------------------------+
251 | oo**(-1+I) | 0 | If the real part of e is negative, then the |
252 | -oo**(-1+I) | | limit is 0. |
253 +--------------+---------+-----------------------------------------------+
255 Because symbolic computations are more flexible than floating point
256 calculations and we prefer to never return an incorrect answer,
257 we choose not to conform to all IEEE 754 conventions. This helps
258 us avoid extra test-case code in the calculation of limits.
260 See Also
261 ========
263 sympy.core.numbers.Infinity
264 sympy.core.numbers.NegativeInfinity
265 sympy.core.numbers.NaN
267 References
268 ==========
270 .. [1] https://en.wikipedia.org/wiki/Exponentiation
271 .. [2] https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero
272 .. [3] https://en.wikipedia.org/wiki/Indeterminate_forms
274 """
275 is_Pow = True
277 __slots__ = ('is_commutative',)
279 args: tuple[Expr, Expr]
280 _args: tuple[Expr, Expr]
282 @cacheit
283 def __new__(cls, b, e, evaluate=None):
284 if evaluate is None:
285 evaluate = global_parameters.evaluate
287 b = _sympify(b)
288 e = _sympify(e)
290 # XXX: This can be removed when non-Expr args are disallowed rather
291 # than deprecated.
292 from .relational import Relational
293 if isinstance(b, Relational) or isinstance(e, Relational):
294 raise TypeError('Relational cannot be used in Pow')
296 # XXX: This should raise TypeError once deprecation period is over:
297 for arg in [b, e]:
298 if not isinstance(arg, Expr):
299 sympy_deprecation_warning(
300 f"""
301 Using non-Expr arguments in Pow is deprecated (in this case, one of the
302 arguments is of type {type(arg).__name__!r}).
304 If you really did intend to construct a power with this base, use the **
305 operator instead.""",
306 deprecated_since_version="1.7",
307 active_deprecations_target="non-expr-args-deprecated",
308 stacklevel=4,
309 )
311 if evaluate:
312 if e is S.ComplexInfinity:
313 return S.NaN
314 if e is S.Infinity:
315 if is_gt(b, S.One):
316 return S.Infinity
317 if is_gt(b, S.NegativeOne) and is_lt(b, S.One):
318 return S.Zero
319 if is_lt(b, S.NegativeOne):
320 if b.is_finite:
321 return S.ComplexInfinity
322 if b.is_finite is False:
323 return S.NaN
324 if e is S.Zero:
325 return S.One
326 elif e is S.One:
327 return b
328 elif e == -1 and not b:
329 return S.ComplexInfinity
330 elif e.__class__.__name__ == "AccumulationBounds":
331 if b == S.Exp1:
332 from sympy.calculus.accumulationbounds import AccumBounds
333 return AccumBounds(Pow(b, e.min), Pow(b, e.max))
334 # autosimplification if base is a number and exp odd/even
335 # if base is Number then the base will end up positive; we
336 # do not do this with arbitrary expressions since symbolic
337 # cancellation might occur as in (x - 1)/(1 - x) -> -1. If
338 # we returned Piecewise((-1, Ne(x, 1))) for such cases then
339 # we could do this...but we don't
340 elif (e.is_Symbol and e.is_integer or e.is_Integer
341 ) and (b.is_number and b.is_Mul or b.is_Number
342 ) and b.could_extract_minus_sign():
343 if e.is_even:
344 b = -b
345 elif e.is_odd:
346 return -Pow(-b, e)
347 if S.NaN in (b, e): # XXX S.NaN**x -> S.NaN under assumption that x != 0
348 return S.NaN
349 elif b is S.One:
350 if abs(e).is_infinite:
351 return S.NaN
352 return S.One
353 else:
354 # recognize base as E
355 from sympy.functions.elementary.exponential import exp_polar
356 if not e.is_Atom and b is not S.Exp1 and not isinstance(b, exp_polar):
357 from .exprtools import factor_terms
358 from sympy.functions.elementary.exponential import log
359 from sympy.simplify.radsimp import fraction
360 c, ex = factor_terms(e, sign=False).as_coeff_Mul()
361 num, den = fraction(ex)
362 if isinstance(den, log) and den.args[0] == b:
363 return S.Exp1**(c*num)
364 elif den.is_Add:
365 from sympy.functions.elementary.complexes import sign, im
366 s = sign(im(b))
367 if s.is_Number and s and den == \
368 log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
369 return S.Exp1**(c*num)
371 obj = b._eval_power(e)
372 if obj is not None:
373 return obj
374 obj = Expr.__new__(cls, b, e)
375 obj = cls._exec_constructor_postprocessors(obj)
376 if not isinstance(obj, Pow):
377 return obj
378 obj.is_commutative = (b.is_commutative and e.is_commutative)
379 return obj
381 def inverse(self, argindex=1):
382 if self.base == S.Exp1:
383 from sympy.functions.elementary.exponential import log
384 return log
385 return None
387 @property
388 def base(self) -> Expr:
389 return self._args[0]
391 @property
392 def exp(self) -> Expr:
393 return self._args[1]
395 @property
396 def kind(self):
397 if self.exp.kind is NumberKind:
398 return self.base.kind
399 else:
400 return UndefinedKind
402 @classmethod
403 def class_key(cls):
404 return 3, 2, cls.__name__
406 def _eval_refine(self, assumptions):
407 from sympy.assumptions.ask import ask, Q
408 b, e = self.as_base_exp()
409 if ask(Q.integer(e), assumptions) and b.could_extract_minus_sign():
410 if ask(Q.even(e), assumptions):
411 return Pow(-b, e)
412 elif ask(Q.odd(e), assumptions):
413 return -Pow(-b, e)
415 def _eval_power(self, other):
416 b, e = self.as_base_exp()
417 if b is S.NaN:
418 return (b**e)**other # let __new__ handle it
420 s = None
421 if other.is_integer:
422 s = 1
423 elif b.is_polar: # e.g. exp_polar, besselj, var('p', polar=True)...
424 s = 1
425 elif e.is_extended_real is not None:
426 from sympy.functions.elementary.complexes import arg, im, re, sign
427 from sympy.functions.elementary.exponential import exp, log
428 from sympy.functions.elementary.integers import floor
429 # helper functions ===========================
430 def _half(e):
431 """Return True if the exponent has a literal 2 as the
432 denominator, else None."""
433 if getattr(e, 'q', None) == 2:
434 return True
435 n, d = e.as_numer_denom()
436 if n.is_integer and d == 2:
437 return True
438 def _n2(e):
439 """Return ``e`` evaluated to a Number with 2 significant
440 digits, else None."""
441 try:
442 rv = e.evalf(2, strict=True)
443 if rv.is_Number:
444 return rv
445 except PrecisionExhausted:
446 pass
447 # ===================================================
448 if e.is_extended_real:
449 # we need _half(other) with constant floor or
450 # floor(S.Half - e*arg(b)/2/pi) == 0
453 # handle -1 as special case
454 if e == -1:
455 # floor arg. is 1/2 + arg(b)/2/pi
456 if _half(other):
457 if b.is_negative is True:
458 return S.NegativeOne**other*Pow(-b, e*other)
459 elif b.is_negative is False: # XXX ok if im(b) != 0?
460 return Pow(b, -other)
461 elif e.is_even:
462 if b.is_extended_real:
463 b = abs(b)
464 if b.is_imaginary:
465 b = abs(im(b))*S.ImaginaryUnit
467 if (abs(e) < 1) == True or e == 1:
468 s = 1 # floor = 0
469 elif b.is_extended_nonnegative:
470 s = 1 # floor = 0
471 elif re(b).is_extended_nonnegative and (abs(e) < 2) == True:
472 s = 1 # floor = 0
473 elif _half(other):
474 s = exp(2*S.Pi*S.ImaginaryUnit*other*floor(
475 S.Half - e*arg(b)/(2*S.Pi)))
476 if s.is_extended_real and _n2(sign(s) - s) == 0:
477 s = sign(s)
478 else:
479 s = None
480 else:
481 # e.is_extended_real is False requires:
482 # _half(other) with constant floor or
483 # floor(S.Half - im(e*log(b))/2/pi) == 0
484 try:
485 s = exp(2*S.ImaginaryUnit*S.Pi*other*
486 floor(S.Half - im(e*log(b))/2/S.Pi))
487 # be careful to test that s is -1 or 1 b/c sign(I) == I:
488 # so check that s is real
489 if s.is_extended_real and _n2(sign(s) - s) == 0:
490 s = sign(s)
491 else:
492 s = None
493 except PrecisionExhausted:
494 s = None
496 if s is not None:
497 return s*Pow(b, e*other)
499 def _eval_Mod(self, q):
500 r"""A dispatched function to compute `b^e \bmod q`, dispatched
501 by ``Mod``.
503 Notes
504 =====
506 Algorithms:
508 1. For unevaluated integer power, use built-in ``pow`` function
509 with 3 arguments, if powers are not too large wrt base.
511 2. For very large powers, use totient reduction if $e \ge \log(m)$.
512 Bound on m, is for safe factorization memory wise i.e. $m^{1/4}$.
513 For pollard-rho to be faster than built-in pow $\log(e) > m^{1/4}$
514 check is added.
516 3. For any unevaluated power found in `b` or `e`, the step 2
517 will be recursed down to the base and the exponent
518 such that the $b \bmod q$ becomes the new base and
519 $\phi(q) + e \bmod \phi(q)$ becomes the new exponent, and then
520 the computation for the reduced expression can be done.
521 """
523 base, exp = self.base, self.exp
525 if exp.is_integer and exp.is_positive:
526 if q.is_integer and base % q == 0:
527 return S.Zero
529 from sympy.ntheory.factor_ import totient
531 if base.is_Integer and exp.is_Integer and q.is_Integer:
532 b, e, m = int(base), int(exp), int(q)
533 mb = m.bit_length()
534 if mb <= 80 and e >= mb and e.bit_length()**4 >= m:
535 phi = int(totient(m))
536 return Integer(pow(b, phi + e%phi, m))
537 return Integer(pow(b, e, m))
539 from .mod import Mod
541 if isinstance(base, Pow) and base.is_integer and base.is_number:
542 base = Mod(base, q)
543 return Mod(Pow(base, exp, evaluate=False), q)
545 if isinstance(exp, Pow) and exp.is_integer and exp.is_number:
546 bit_length = int(q).bit_length()
547 # XXX Mod-Pow actually attempts to do a hanging evaluation
548 # if this dispatched function returns None.
549 # May need some fixes in the dispatcher itself.
550 if bit_length <= 80:
551 phi = totient(q)
552 exp = phi + Mod(exp, phi)
553 return Mod(Pow(base, exp, evaluate=False), q)
555 def _eval_is_even(self):
556 if self.exp.is_integer and self.exp.is_positive:
557 return self.base.is_even
559 def _eval_is_negative(self):
560 ext_neg = Pow._eval_is_extended_negative(self)
561 if ext_neg is True:
562 return self.is_finite
563 return ext_neg
565 def _eval_is_extended_positive(self):
566 if self.base == self.exp:
567 if self.base.is_extended_nonnegative:
568 return True
569 elif self.base.is_positive:
570 if self.exp.is_real:
571 return True
572 elif self.base.is_extended_negative:
573 if self.exp.is_even:
574 return True
575 if self.exp.is_odd:
576 return False
577 elif self.base.is_zero:
578 if self.exp.is_extended_real:
579 return self.exp.is_zero
580 elif self.base.is_extended_nonpositive:
581 if self.exp.is_odd:
582 return False
583 elif self.base.is_imaginary:
584 if self.exp.is_integer:
585 m = self.exp % 4
586 if m.is_zero:
587 return True
588 if m.is_integer and m.is_zero is False:
589 return False
590 if self.exp.is_imaginary:
591 from sympy.functions.elementary.exponential import log
592 return log(self.base).is_imaginary
594 def _eval_is_extended_negative(self):
595 if self.exp is S.Half:
596 if self.base.is_complex or self.base.is_extended_real:
597 return False
598 if self.base.is_extended_negative:
599 if self.exp.is_odd and self.base.is_finite:
600 return True
601 if self.exp.is_even:
602 return False
603 elif self.base.is_extended_positive:
604 if self.exp.is_extended_real:
605 return False
606 elif self.base.is_zero:
607 if self.exp.is_extended_real:
608 return False
609 elif self.base.is_extended_nonnegative:
610 if self.exp.is_extended_nonnegative:
611 return False
612 elif self.base.is_extended_nonpositive:
613 if self.exp.is_even:
614 return False
615 elif self.base.is_extended_real:
616 if self.exp.is_even:
617 return False
619 def _eval_is_zero(self):
620 if self.base.is_zero:
621 if self.exp.is_extended_positive:
622 return True
623 elif self.exp.is_extended_nonpositive:
624 return False
625 elif self.base == S.Exp1:
626 return self.exp is S.NegativeInfinity
627 elif self.base.is_zero is False:
628 if self.base.is_finite and self.exp.is_finite:
629 return False
630 elif self.exp.is_negative:
631 return self.base.is_infinite
632 elif self.exp.is_nonnegative:
633 return False
634 elif self.exp.is_infinite and self.exp.is_extended_real:
635 if (1 - abs(self.base)).is_extended_positive:
636 return self.exp.is_extended_positive
637 elif (1 - abs(self.base)).is_extended_negative:
638 return self.exp.is_extended_negative
639 elif self.base.is_finite and self.exp.is_negative:
640 # when self.base.is_zero is None
641 return False
643 def _eval_is_integer(self):
644 b, e = self.args
645 if b.is_rational:
646 if b.is_integer is False and e.is_positive:
647 return False # rat**nonneg
648 if b.is_integer and e.is_integer:
649 if b is S.NegativeOne:
650 return True
651 if e.is_nonnegative or e.is_positive:
652 return True
653 if b.is_integer and e.is_negative and (e.is_finite or e.is_integer):
654 if fuzzy_not((b - 1).is_zero) and fuzzy_not((b + 1).is_zero):
655 return False
656 if b.is_Number and e.is_Number:
657 check = self.func(*self.args)
658 return check.is_Integer
659 if e.is_negative and b.is_positive and (b - 1).is_positive:
660 return False
661 if e.is_negative and b.is_negative and (b + 1).is_negative:
662 return False
664 def _eval_is_extended_real(self):
665 if self.base is S.Exp1:
666 if self.exp.is_extended_real:
667 return True
668 elif self.exp.is_imaginary:
669 return (2*S.ImaginaryUnit*self.exp/S.Pi).is_even
671 from sympy.functions.elementary.exponential import log, exp
672 real_b = self.base.is_extended_real
673 if real_b is None:
674 if self.base.func == exp and self.base.exp.is_imaginary:
675 return self.exp.is_imaginary
676 if self.base.func == Pow and self.base.base is S.Exp1 and self.base.exp.is_imaginary:
677 return self.exp.is_imaginary
678 return
679 real_e = self.exp.is_extended_real
680 if real_e is None:
681 return
682 if real_b and real_e:
683 if self.base.is_extended_positive:
684 return True
685 elif self.base.is_extended_nonnegative and self.exp.is_extended_nonnegative:
686 return True
687 elif self.exp.is_integer and self.base.is_extended_nonzero:
688 return True
689 elif self.exp.is_integer and self.exp.is_nonnegative:
690 return True
691 elif self.base.is_extended_negative:
692 if self.exp.is_Rational:
693 return False
694 if real_e and self.exp.is_extended_negative and self.base.is_zero is False:
695 return Pow(self.base, -self.exp).is_extended_real
696 im_b = self.base.is_imaginary
697 im_e = self.exp.is_imaginary
698 if im_b:
699 if self.exp.is_integer:
700 if self.exp.is_even:
701 return True
702 elif self.exp.is_odd:
703 return False
704 elif im_e and log(self.base).is_imaginary:
705 return True
706 elif self.exp.is_Add:
707 c, a = self.exp.as_coeff_Add()
708 if c and c.is_Integer:
709 return Mul(
710 self.base**c, self.base**a, evaluate=False).is_extended_real
711 elif self.base in (-S.ImaginaryUnit, S.ImaginaryUnit):
712 if (self.exp/2).is_integer is False:
713 return False
714 if real_b and im_e:
715 if self.base is S.NegativeOne:
716 return True
717 c = self.exp.coeff(S.ImaginaryUnit)
718 if c:
719 if self.base.is_rational and c.is_rational:
720 if self.base.is_nonzero and (self.base - 1).is_nonzero and c.is_nonzero:
721 return False
722 ok = (c*log(self.base)/S.Pi).is_integer
723 if ok is not None:
724 return ok
726 if real_b is False and real_e: # we already know it's not imag
727 from sympy.functions.elementary.complexes import arg
728 i = arg(self.base)*self.exp/S.Pi
729 if i.is_complex: # finite
730 return i.is_integer
732 def _eval_is_complex(self):
734 if self.base == S.Exp1:
735 return fuzzy_or([self.exp.is_complex, self.exp.is_extended_negative])
737 if all(a.is_complex for a in self.args) and self._eval_is_finite():
738 return True
740 def _eval_is_imaginary(self):
741 if self.base.is_commutative is False:
742 return False
744 if self.base.is_imaginary:
745 if self.exp.is_integer:
746 odd = self.exp.is_odd
747 if odd is not None:
748 return odd
749 return
751 if self.base == S.Exp1:
752 f = 2 * self.exp / (S.Pi*S.ImaginaryUnit)
753 # exp(pi*integer) = 1 or -1, so not imaginary
754 if f.is_even:
755 return False
756 # exp(pi*integer + pi/2) = I or -I, so it is imaginary
757 if f.is_odd:
758 return True
759 return None
761 if self.exp.is_imaginary:
762 from sympy.functions.elementary.exponential import log
763 imlog = log(self.base).is_imaginary
764 if imlog is not None:
765 return False # I**i -> real; (2*I)**i -> complex ==> not imaginary
767 if self.base.is_extended_real and self.exp.is_extended_real:
768 if self.base.is_positive:
769 return False
770 else:
771 rat = self.exp.is_rational
772 if not rat:
773 return rat
774 if self.exp.is_integer:
775 return False
776 else:
777 half = (2*self.exp).is_integer
778 if half:
779 return self.base.is_negative
780 return half
782 if self.base.is_extended_real is False: # we already know it's not imag
783 from sympy.functions.elementary.complexes import arg
784 i = arg(self.base)*self.exp/S.Pi
785 isodd = (2*i).is_odd
786 if isodd is not None:
787 return isodd
789 def _eval_is_odd(self):
790 if self.exp.is_integer:
791 if self.exp.is_positive:
792 return self.base.is_odd
793 elif self.exp.is_nonnegative and self.base.is_odd:
794 return True
795 elif self.base is S.NegativeOne:
796 return True
798 def _eval_is_finite(self):
799 if self.exp.is_negative:
800 if self.base.is_zero:
801 return False
802 if self.base.is_infinite or self.base.is_nonzero:
803 return True
804 c1 = self.base.is_finite
805 if c1 is None:
806 return
807 c2 = self.exp.is_finite
808 if c2 is None:
809 return
810 if c1 and c2:
811 if self.exp.is_nonnegative or fuzzy_not(self.base.is_zero):
812 return True
814 def _eval_is_prime(self):
815 '''
816 An integer raised to the n(>=2)-th power cannot be a prime.
817 '''
818 if self.base.is_integer and self.exp.is_integer and (self.exp - 1).is_positive:
819 return False
821 def _eval_is_composite(self):
822 """
823 A power is composite if both base and exponent are greater than 1
824 """
825 if (self.base.is_integer and self.exp.is_integer and
826 ((self.base - 1).is_positive and (self.exp - 1).is_positive or
827 (self.base + 1).is_negative and self.exp.is_positive and self.exp.is_even)):
828 return True
830 def _eval_is_polar(self):
831 return self.base.is_polar
833 def _eval_subs(self, old, new):
834 from sympy.calculus.accumulationbounds import AccumBounds
836 if isinstance(self.exp, AccumBounds):
837 b = self.base.subs(old, new)
838 e = self.exp.subs(old, new)
839 if isinstance(e, AccumBounds):
840 return e.__rpow__(b)
841 return self.func(b, e)
843 from sympy.functions.elementary.exponential import exp, log
845 def _check(ct1, ct2, old):
846 """Return (bool, pow, remainder_pow) where, if bool is True, then the
847 exponent of Pow `old` will combine with `pow` so the substitution
848 is valid, otherwise bool will be False.
850 For noncommutative objects, `pow` will be an integer, and a factor
851 `Pow(old.base, remainder_pow)` needs to be included. If there is
852 no such factor, None is returned. For commutative objects,
853 remainder_pow is always None.
855 cti are the coefficient and terms of an exponent of self or old
856 In this _eval_subs routine a change like (b**(2*x)).subs(b**x, y)
857 will give y**2 since (b**x)**2 == b**(2*x); if that equality does
858 not hold then the substitution should not occur so `bool` will be
859 False.
861 """
862 coeff1, terms1 = ct1
863 coeff2, terms2 = ct2
864 if terms1 == terms2:
865 if old.is_commutative:
866 # Allow fractional powers for commutative objects
867 pow = coeff1/coeff2
868 try:
869 as_int(pow, strict=False)
870 combines = True
871 except ValueError:
872 b, e = old.as_base_exp()
873 # These conditions ensure that (b**e)**f == b**(e*f) for any f
874 combines = b.is_positive and e.is_real or b.is_nonnegative and e.is_nonnegative
876 return combines, pow, None
877 else:
878 # With noncommutative symbols, substitute only integer powers
879 if not isinstance(terms1, tuple):
880 terms1 = (terms1,)
881 if not all(term.is_integer for term in terms1):
882 return False, None, None
884 try:
885 # Round pow toward zero
886 pow, remainder = divmod(as_int(coeff1), as_int(coeff2))
887 if pow < 0 and remainder != 0:
888 pow += 1
889 remainder -= as_int(coeff2)
891 if remainder == 0:
892 remainder_pow = None
893 else:
894 remainder_pow = Mul(remainder, *terms1)
896 return True, pow, remainder_pow
897 except ValueError:
898 # Can't substitute
899 pass
901 return False, None, None
903 if old == self.base or (old == exp and self.base == S.Exp1):
904 if new.is_Function and isinstance(new, Callable):
905 return new(self.exp._subs(old, new))
906 else:
907 return new**self.exp._subs(old, new)
909 # issue 10829: (4**x - 3*y + 2).subs(2**x, y) -> y**2 - 3*y + 2
910 if isinstance(old, self.func) and self.exp == old.exp:
911 l = log(self.base, old.base)
912 if l.is_Number:
913 return Pow(new, l)
915 if isinstance(old, self.func) and self.base == old.base:
916 if self.exp.is_Add is False:
917 ct1 = self.exp.as_independent(Symbol, as_Add=False)
918 ct2 = old.exp.as_independent(Symbol, as_Add=False)
919 ok, pow, remainder_pow = _check(ct1, ct2, old)
920 if ok:
921 # issue 5180: (x**(6*y)).subs(x**(3*y),z)->z**2
922 result = self.func(new, pow)
923 if remainder_pow is not None:
924 result = Mul(result, Pow(old.base, remainder_pow))
925 return result
926 else: # b**(6*x + a).subs(b**(3*x), y) -> y**2 * b**a
927 # exp(exp(x) + exp(x**2)).subs(exp(exp(x)), w) -> w * exp(exp(x**2))
928 oarg = old.exp
929 new_l = []
930 o_al = []
931 ct2 = oarg.as_coeff_mul()
932 for a in self.exp.args:
933 newa = a._subs(old, new)
934 ct1 = newa.as_coeff_mul()
935 ok, pow, remainder_pow = _check(ct1, ct2, old)
936 if ok:
937 new_l.append(new**pow)
938 if remainder_pow is not None:
939 o_al.append(remainder_pow)
940 continue
941 elif not old.is_commutative and not newa.is_integer:
942 # If any term in the exponent is non-integer,
943 # we do not do any substitutions in the noncommutative case
944 return
945 o_al.append(newa)
946 if new_l:
947 expo = Add(*o_al)
948 new_l.append(Pow(self.base, expo, evaluate=False) if expo != 1 else self.base)
949 return Mul(*new_l)
951 if (isinstance(old, exp) or (old.is_Pow and old.base is S.Exp1)) and self.exp.is_extended_real and self.base.is_positive:
952 ct1 = old.exp.as_independent(Symbol, as_Add=False)
953 ct2 = (self.exp*log(self.base)).as_independent(
954 Symbol, as_Add=False)
955 ok, pow, remainder_pow = _check(ct1, ct2, old)
956 if ok:
957 result = self.func(new, pow) # (2**x).subs(exp(x*log(2)), z) -> z
958 if remainder_pow is not None:
959 result = Mul(result, Pow(old.base, remainder_pow))
960 return result
962 def as_base_exp(self):
963 """Return base and exp of self.
965 Explanation
966 ===========
968 If base a Rational less than 1, then return 1/Rational, -exp.
969 If this extra processing is not needed, the base and exp
970 properties will give the raw arguments.
972 Examples
973 ========
975 >>> from sympy import Pow, S
976 >>> p = Pow(S.Half, 2, evaluate=False)
977 >>> p.as_base_exp()
978 (2, -2)
979 >>> p.args
980 (1/2, 2)
981 >>> p.base, p.exp
982 (1/2, 2)
984 """
986 b, e = self.args
987 if b.is_Rational and b.p < b.q and b.p > 0:
988 return 1/b, -e
989 return b, e
991 def _eval_adjoint(self):
992 from sympy.functions.elementary.complexes import adjoint
993 i, p = self.exp.is_integer, self.base.is_positive
994 if i:
995 return adjoint(self.base)**self.exp
996 if p:
997 return self.base**adjoint(self.exp)
998 if i is False and p is False:
999 expanded = expand_complex(self)
1000 if expanded != self:
1001 return adjoint(expanded)
1003 def _eval_conjugate(self):
1004 from sympy.functions.elementary.complexes import conjugate as c
1005 i, p = self.exp.is_integer, self.base.is_positive
1006 if i:
1007 return c(self.base)**self.exp
1008 if p:
1009 return self.base**c(self.exp)
1010 if i is False and p is False:
1011 expanded = expand_complex(self)
1012 if expanded != self:
1013 return c(expanded)
1014 if self.is_extended_real:
1015 return self
1017 def _eval_transpose(self):
1018 from sympy.functions.elementary.complexes import transpose
1019 if self.base == S.Exp1:
1020 return self.func(S.Exp1, self.exp.transpose())
1021 i, p = self.exp.is_integer, (self.base.is_complex or self.base.is_infinite)
1022 if p:
1023 return self.base**self.exp
1024 if i:
1025 return transpose(self.base)**self.exp
1026 if i is False and p is False:
1027 expanded = expand_complex(self)
1028 if expanded != self:
1029 return transpose(expanded)
1031 def _eval_expand_power_exp(self, **hints):
1032 """a**(n + m) -> a**n*a**m"""
1033 b = self.base
1034 e = self.exp
1035 if b == S.Exp1:
1036 from sympy.concrete.summations import Sum
1037 if isinstance(e, Sum) and e.is_commutative:
1038 from sympy.concrete.products import Product
1039 return Product(self.func(b, e.function), *e.limits)
1040 if e.is_Add and (hints.get('force', False) or
1041 b.is_zero is False or e._all_nonneg_or_nonppos()):
1042 if e.is_commutative:
1043 return Mul(*[self.func(b, x) for x in e.args])
1044 if b.is_commutative:
1045 c, nc = sift(e.args, lambda x: x.is_commutative, binary=True)
1046 if c:
1047 return Mul(*[self.func(b, x) for x in c]
1048 )*b**Add._from_args(nc)
1049 return self
1051 def _eval_expand_power_base(self, **hints):
1052 """(a*b)**n -> a**n * b**n"""
1053 force = hints.get('force', False)
1055 b = self.base
1056 e = self.exp
1057 if not b.is_Mul:
1058 return self
1060 cargs, nc = b.args_cnc(split_1=False)
1062 # expand each term - this is top-level-only
1063 # expansion but we have to watch out for things
1064 # that don't have an _eval_expand method
1065 if nc:
1066 nc = [i._eval_expand_power_base(**hints)
1067 if hasattr(i, '_eval_expand_power_base') else i
1068 for i in nc]
1070 if e.is_Integer:
1071 if e.is_positive:
1072 rv = Mul(*nc*e)
1073 else:
1074 rv = Mul(*[i**-1 for i in nc[::-1]]*-e)
1075 if cargs:
1076 rv *= Mul(*cargs)**e
1077 return rv
1079 if not cargs:
1080 return self.func(Mul(*nc), e, evaluate=False)
1082 nc = [Mul(*nc)]
1084 # sift the commutative bases
1085 other, maybe_real = sift(cargs, lambda x: x.is_extended_real is False,
1086 binary=True)
1087 def pred(x):
1088 if x is S.ImaginaryUnit:
1089 return S.ImaginaryUnit
1090 polar = x.is_polar
1091 if polar:
1092 return True
1093 if polar is None:
1094 return fuzzy_bool(x.is_extended_nonnegative)
1095 sifted = sift(maybe_real, pred)
1096 nonneg = sifted[True]
1097 other += sifted[None]
1098 neg = sifted[False]
1099 imag = sifted[S.ImaginaryUnit]
1100 if imag:
1101 I = S.ImaginaryUnit
1102 i = len(imag) % 4
1103 if i == 0:
1104 pass
1105 elif i == 1:
1106 other.append(I)
1107 elif i == 2:
1108 if neg:
1109 nonn = -neg.pop()
1110 if nonn is not S.One:
1111 nonneg.append(nonn)
1112 else:
1113 neg.append(S.NegativeOne)
1114 else:
1115 if neg:
1116 nonn = -neg.pop()
1117 if nonn is not S.One:
1118 nonneg.append(nonn)
1119 else:
1120 neg.append(S.NegativeOne)
1121 other.append(I)
1122 del imag
1124 # bring out the bases that can be separated from the base
1126 if force or e.is_integer:
1127 # treat all commutatives the same and put nc in other
1128 cargs = nonneg + neg + other
1129 other = nc
1130 else:
1131 # this is just like what is happening automatically, except
1132 # that now we are doing it for an arbitrary exponent for which
1133 # no automatic expansion is done
1135 assert not e.is_Integer
1137 # handle negatives by making them all positive and putting
1138 # the residual -1 in other
1139 if len(neg) > 1:
1140 o = S.One
1141 if not other and neg[0].is_Number:
1142 o *= neg.pop(0)
1143 if len(neg) % 2:
1144 o = -o
1145 for n in neg:
1146 nonneg.append(-n)
1147 if o is not S.One:
1148 other.append(o)
1149 elif neg and other:
1150 if neg[0].is_Number and neg[0] is not S.NegativeOne:
1151 other.append(S.NegativeOne)
1152 nonneg.append(-neg[0])
1153 else:
1154 other.extend(neg)
1155 else:
1156 other.extend(neg)
1157 del neg
1159 cargs = nonneg
1160 other += nc
1162 rv = S.One
1163 if cargs:
1164 if e.is_Rational:
1165 npow, cargs = sift(cargs, lambda x: x.is_Pow and
1166 x.exp.is_Rational and x.base.is_number,
1167 binary=True)
1168 rv = Mul(*[self.func(b.func(*b.args), e) for b in npow])
1169 rv *= Mul(*[self.func(b, e, evaluate=False) for b in cargs])
1170 if other:
1171 rv *= self.func(Mul(*other), e, evaluate=False)
1172 return rv
1174 def _eval_expand_multinomial(self, **hints):
1175 """(a + b + ..)**n -> a**n + n*a**(n-1)*b + .., n is nonzero integer"""
1177 base, exp = self.args
1178 result = self
1180 if exp.is_Rational and exp.p > 0 and base.is_Add:
1181 if not exp.is_Integer:
1182 n = Integer(exp.p // exp.q)
1184 if not n:
1185 return result
1186 else:
1187 radical, result = self.func(base, exp - n), []
1189 expanded_base_n = self.func(base, n)
1190 if expanded_base_n.is_Pow:
1191 expanded_base_n = \
1192 expanded_base_n._eval_expand_multinomial()
1193 for term in Add.make_args(expanded_base_n):
1194 result.append(term*radical)
1196 return Add(*result)
1198 n = int(exp)
1200 if base.is_commutative:
1201 order_terms, other_terms = [], []
1203 for b in base.args:
1204 if b.is_Order:
1205 order_terms.append(b)
1206 else:
1207 other_terms.append(b)
1209 if order_terms:
1210 # (f(x) + O(x^n))^m -> f(x)^m + m*f(x)^{m-1} *O(x^n)
1211 f = Add(*other_terms)
1212 o = Add(*order_terms)
1214 if n == 2:
1215 return expand_multinomial(f**n, deep=False) + n*f*o
1216 else:
1217 g = expand_multinomial(f**(n - 1), deep=False)
1218 return expand_mul(f*g, deep=False) + n*g*o
1220 if base.is_number:
1221 # Efficiently expand expressions of the form (a + b*I)**n
1222 # where 'a' and 'b' are real numbers and 'n' is integer.
1223 a, b = base.as_real_imag()
1225 if a.is_Rational and b.is_Rational:
1226 if not a.is_Integer:
1227 if not b.is_Integer:
1228 k = self.func(a.q * b.q, n)
1229 a, b = a.p*b.q, a.q*b.p
1230 else:
1231 k = self.func(a.q, n)
1232 a, b = a.p, a.q*b
1233 elif not b.is_Integer:
1234 k = self.func(b.q, n)
1235 a, b = a*b.q, b.p
1236 else:
1237 k = 1
1239 a, b, c, d = int(a), int(b), 1, 0
1241 while n:
1242 if n & 1:
1243 c, d = a*c - b*d, b*c + a*d
1244 n -= 1
1245 a, b = a*a - b*b, 2*a*b
1246 n //= 2
1248 I = S.ImaginaryUnit
1250 if k == 1:
1251 return c + I*d
1252 else:
1253 return Integer(c)/k + I*d/k
1255 p = other_terms
1256 # (x + y)**3 -> x**3 + 3*x**2*y + 3*x*y**2 + y**3
1257 # in this particular example:
1258 # p = [x,y]; n = 3
1259 # so now it's easy to get the correct result -- we get the
1260 # coefficients first:
1261 from sympy.ntheory.multinomial import multinomial_coefficients
1262 from sympy.polys.polyutils import basic_from_dict
1263 expansion_dict = multinomial_coefficients(len(p), n)
1264 # in our example: {(3, 0): 1, (1, 2): 3, (0, 3): 1, (2, 1): 3}
1265 # and now construct the expression.
1266 return basic_from_dict(expansion_dict, *p)
1267 else:
1268 if n == 2:
1269 return Add(*[f*g for f in base.args for g in base.args])
1270 else:
1271 multi = (base**(n - 1))._eval_expand_multinomial()
1272 if multi.is_Add:
1273 return Add(*[f*g for f in base.args
1274 for g in multi.args])
1275 else:
1276 # XXX can this ever happen if base was an Add?
1277 return Add(*[f*multi for f in base.args])
1278 elif (exp.is_Rational and exp.p < 0 and base.is_Add and
1279 abs(exp.p) > exp.q):
1280 return 1 / self.func(base, -exp)._eval_expand_multinomial()
1281 elif exp.is_Add and base.is_Number and (hints.get('force', False) or
1282 base.is_zero is False or exp._all_nonneg_or_nonppos()):
1283 # a + b a b
1284 # n --> n n, where n, a, b are Numbers
1285 # XXX should be in expand_power_exp?
1286 coeff, tail = [], []
1287 for term in exp.args:
1288 if term.is_Number:
1289 coeff.append(self.func(base, term))
1290 else:
1291 tail.append(term)
1292 return Mul(*(coeff + [self.func(base, Add._from_args(tail))]))
1293 else:
1294 return result
1296 def as_real_imag(self, deep=True, **hints):
1297 if self.exp.is_Integer:
1298 from sympy.polys.polytools import poly
1300 exp = self.exp
1301 re_e, im_e = self.base.as_real_imag(deep=deep)
1302 if not im_e:
1303 return self, S.Zero
1304 a, b = symbols('a b', cls=Dummy)
1305 if exp >= 0:
1306 if re_e.is_Number and im_e.is_Number:
1307 # We can be more efficient in this case
1308 expr = expand_multinomial(self.base**exp)
1309 if expr != self:
1310 return expr.as_real_imag()
1312 expr = poly(
1313 (a + b)**exp) # a = re, b = im; expr = (a + b*I)**exp
1314 else:
1315 mag = re_e**2 + im_e**2
1316 re_e, im_e = re_e/mag, -im_e/mag
1317 if re_e.is_Number and im_e.is_Number:
1318 # We can be more efficient in this case
1319 expr = expand_multinomial((re_e + im_e*S.ImaginaryUnit)**-exp)
1320 if expr != self:
1321 return expr.as_real_imag()
1323 expr = poly((a + b)**-exp)
1325 # Terms with even b powers will be real
1326 r = [i for i in expr.terms() if not i[0][1] % 2]
1327 re_part = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r])
1328 # Terms with odd b powers will be imaginary
1329 r = [i for i in expr.terms() if i[0][1] % 4 == 1]
1330 im_part1 = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r])
1331 r = [i for i in expr.terms() if i[0][1] % 4 == 3]
1332 im_part3 = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r])
1334 return (re_part.subs({a: re_e, b: S.ImaginaryUnit*im_e}),
1335 im_part1.subs({a: re_e, b: im_e}) + im_part3.subs({a: re_e, b: -im_e}))
1337 from sympy.functions.elementary.trigonometric import atan2, cos, sin
1339 if self.exp.is_Rational:
1340 re_e, im_e = self.base.as_real_imag(deep=deep)
1342 if im_e.is_zero and self.exp is S.Half:
1343 if re_e.is_extended_nonnegative:
1344 return self, S.Zero
1345 if re_e.is_extended_nonpositive:
1346 return S.Zero, (-self.base)**self.exp
1348 # XXX: This is not totally correct since for x**(p/q) with
1349 # x being imaginary there are actually q roots, but
1350 # only a single one is returned from here.
1351 r = self.func(self.func(re_e, 2) + self.func(im_e, 2), S.Half)
1353 t = atan2(im_e, re_e)
1355 rp, tp = self.func(r, self.exp), t*self.exp
1357 return rp*cos(tp), rp*sin(tp)
1358 elif self.base is S.Exp1:
1359 from sympy.functions.elementary.exponential import exp
1360 re_e, im_e = self.exp.as_real_imag()
1361 if deep:
1362 re_e = re_e.expand(deep, **hints)
1363 im_e = im_e.expand(deep, **hints)
1364 c, s = cos(im_e), sin(im_e)
1365 return exp(re_e)*c, exp(re_e)*s
1366 else:
1367 from sympy.functions.elementary.complexes import im, re
1368 if deep:
1369 hints['complex'] = False
1371 expanded = self.expand(deep, **hints)
1372 if hints.get('ignore') == expanded:
1373 return None
1374 else:
1375 return (re(expanded), im(expanded))
1376 else:
1377 return re(self), im(self)
1379 def _eval_derivative(self, s):
1380 from sympy.functions.elementary.exponential import log
1381 dbase = self.base.diff(s)
1382 dexp = self.exp.diff(s)
1383 return self * (dexp * log(self.base) + dbase * self.exp/self.base)
1385 def _eval_evalf(self, prec):
1386 base, exp = self.as_base_exp()
1387 if base == S.Exp1:
1388 # Use mpmath function associated to class "exp":
1389 from sympy.functions.elementary.exponential import exp as exp_function
1390 return exp_function(self.exp, evaluate=False)._eval_evalf(prec)
1391 base = base._evalf(prec)
1392 if not exp.is_Integer:
1393 exp = exp._evalf(prec)
1394 if exp.is_negative and base.is_number and base.is_extended_real is False:
1395 base = base.conjugate() / (base * base.conjugate())._evalf(prec)
1396 exp = -exp
1397 return self.func(base, exp).expand()
1398 return self.func(base, exp)
1400 def _eval_is_polynomial(self, syms):
1401 if self.exp.has(*syms):
1402 return False
1404 if self.base.has(*syms):
1405 return bool(self.base._eval_is_polynomial(syms) and
1406 self.exp.is_Integer and (self.exp >= 0))
1407 else:
1408 return True
1410 def _eval_is_rational(self):
1411 # The evaluation of self.func below can be very expensive in the case
1412 # of integer**integer if the exponent is large. We should try to exit
1413 # before that if possible:
1414 if (self.exp.is_integer and self.base.is_rational
1415 and fuzzy_not(fuzzy_and([self.exp.is_negative, self.base.is_zero]))):
1416 return True
1417 p = self.func(*self.as_base_exp()) # in case it's unevaluated
1418 if not p.is_Pow:
1419 return p.is_rational
1420 b, e = p.as_base_exp()
1421 if e.is_Rational and b.is_Rational:
1422 # we didn't check that e is not an Integer
1423 # because Rational**Integer autosimplifies
1424 return False
1425 if e.is_integer:
1426 if b.is_rational:
1427 if fuzzy_not(b.is_zero) or e.is_nonnegative:
1428 return True
1429 if b == e: # always rational, even for 0**0
1430 return True
1431 elif b.is_irrational:
1432 return e.is_zero
1433 if b is S.Exp1:
1434 if e.is_rational and e.is_nonzero:
1435 return False
1437 def _eval_is_algebraic(self):
1438 def _is_one(expr):
1439 try:
1440 return (expr - 1).is_zero
1441 except ValueError:
1442 # when the operation is not allowed
1443 return False
1445 if self.base.is_zero or _is_one(self.base):
1446 return True
1447 elif self.base is S.Exp1:
1448 s = self.func(*self.args)
1449 if s.func == self.func:
1450 if self.exp.is_nonzero:
1451 if self.exp.is_algebraic:
1452 return False
1453 elif (self.exp/S.Pi).is_rational:
1454 return False
1455 elif (self.exp/(S.ImaginaryUnit*S.Pi)).is_rational:
1456 return True
1457 else:
1458 return s.is_algebraic
1459 elif self.exp.is_rational:
1460 if self.base.is_algebraic is False:
1461 return self.exp.is_zero
1462 if self.base.is_zero is False:
1463 if self.exp.is_nonzero:
1464 return self.base.is_algebraic
1465 elif self.base.is_algebraic:
1466 return True
1467 if self.exp.is_positive:
1468 return self.base.is_algebraic
1469 elif self.base.is_algebraic and self.exp.is_algebraic:
1470 if ((fuzzy_not(self.base.is_zero)
1471 and fuzzy_not(_is_one(self.base)))
1472 or self.base.is_integer is False
1473 or self.base.is_irrational):
1474 return self.exp.is_rational
1476 def _eval_is_rational_function(self, syms):
1477 if self.exp.has(*syms):
1478 return False
1480 if self.base.has(*syms):
1481 return self.base._eval_is_rational_function(syms) and \
1482 self.exp.is_Integer
1483 else:
1484 return True
1486 def _eval_is_meromorphic(self, x, a):
1487 # f**g is meromorphic if g is an integer and f is meromorphic.
1488 # E**(log(f)*g) is meromorphic if log(f)*g is meromorphic
1489 # and finite.
1490 base_merom = self.base._eval_is_meromorphic(x, a)
1491 exp_integer = self.exp.is_Integer
1492 if exp_integer:
1493 return base_merom
1495 exp_merom = self.exp._eval_is_meromorphic(x, a)
1496 if base_merom is False:
1497 # f**g = E**(log(f)*g) may be meromorphic if the
1498 # singularities of log(f) and g cancel each other,
1499 # for example, if g = 1/log(f). Hence,
1500 return False if exp_merom else None
1501 elif base_merom is None:
1502 return None
1504 b = self.base.subs(x, a)
1505 # b is extended complex as base is meromorphic.
1506 # log(base) is finite and meromorphic when b != 0, zoo.
1507 b_zero = b.is_zero
1508 if b_zero:
1509 log_defined = False
1510 else:
1511 log_defined = fuzzy_and((b.is_finite, fuzzy_not(b_zero)))
1513 if log_defined is False: # zero or pole of base
1514 return exp_integer # False or None
1515 elif log_defined is None:
1516 return None
1518 if not exp_merom:
1519 return exp_merom # False or None
1521 return self.exp.subs(x, a).is_finite
1523 def _eval_is_algebraic_expr(self, syms):
1524 if self.exp.has(*syms):
1525 return False
1527 if self.base.has(*syms):
1528 return self.base._eval_is_algebraic_expr(syms) and \
1529 self.exp.is_Rational
1530 else:
1531 return True
1533 def _eval_rewrite_as_exp(self, base, expo, **kwargs):
1534 from sympy.functions.elementary.exponential import exp, log
1536 if base.is_zero or base.has(exp) or expo.has(exp):
1537 return base**expo
1539 if base.has(Symbol):
1540 # delay evaluation if expo is non symbolic
1541 # (as exp(x*log(5)) automatically reduces to x**5)
1542 if global_parameters.exp_is_pow:
1543 return Pow(S.Exp1, log(base)*expo, evaluate=expo.has(Symbol))
1544 else:
1545 return exp(log(base)*expo, evaluate=expo.has(Symbol))
1547 else:
1548 from sympy.functions.elementary.complexes import arg, Abs
1549 return exp((log(Abs(base)) + S.ImaginaryUnit*arg(base))*expo)
1551 def as_numer_denom(self):
1552 if not self.is_commutative:
1553 return self, S.One
1554 base, exp = self.as_base_exp()
1555 n, d = base.as_numer_denom()
1556 # this should be the same as ExpBase.as_numer_denom wrt
1557 # exponent handling
1558 neg_exp = exp.is_negative
1559 if exp.is_Mul and not neg_exp and not exp.is_positive:
1560 neg_exp = exp.could_extract_minus_sign()
1561 int_exp = exp.is_integer
1562 # the denominator cannot be separated from the numerator if
1563 # its sign is unknown unless the exponent is an integer, e.g.
1564 # sqrt(a/b) != sqrt(a)/sqrt(b) when a=1 and b=-1. But if the
1565 # denominator is negative the numerator and denominator can
1566 # be negated and the denominator (now positive) separated.
1567 if not (d.is_extended_real or int_exp):
1568 n = base
1569 d = S.One
1570 dnonpos = d.is_nonpositive
1571 if dnonpos:
1572 n, d = -n, -d
1573 elif dnonpos is None and not int_exp:
1574 n = base
1575 d = S.One
1576 if neg_exp:
1577 n, d = d, n
1578 exp = -exp
1579 if exp.is_infinite:
1580 if n is S.One and d is not S.One:
1581 return n, self.func(d, exp)
1582 if n is not S.One and d is S.One:
1583 return self.func(n, exp), d
1584 return self.func(n, exp), self.func(d, exp)
1586 def matches(self, expr, repl_dict=None, old=False):
1587 expr = _sympify(expr)
1588 if repl_dict is None:
1589 repl_dict = {}
1591 # special case, pattern = 1 and expr.exp can match to 0
1592 if expr is S.One:
1593 d = self.exp.matches(S.Zero, repl_dict)
1594 if d is not None:
1595 return d
1597 # make sure the expression to be matched is an Expr
1598 if not isinstance(expr, Expr):
1599 return None
1601 b, e = expr.as_base_exp()
1603 # special case number
1604 sb, se = self.as_base_exp()
1605 if sb.is_Symbol and se.is_Integer and expr:
1606 if e.is_rational:
1607 return sb.matches(b**(e/se), repl_dict)
1608 return sb.matches(expr**(1/se), repl_dict)
1610 d = repl_dict.copy()
1611 d = self.base.matches(b, d)
1612 if d is None:
1613 return None
1615 d = self.exp.xreplace(d).matches(e, d)
1616 if d is None:
1617 return Expr.matches(self, expr, repl_dict)
1618 return d
1620 def _eval_nseries(self, x, n, logx, cdir=0):
1621 # NOTE! This function is an important part of the gruntz algorithm
1622 # for computing limits. It has to return a generalized power
1623 # series with coefficients in C(log, log(x)). In more detail:
1624 # It has to return an expression
1625 # c_0*x**e_0 + c_1*x**e_1 + ... (finitely many terms)
1626 # where e_i are numbers (not necessarily integers) and c_i are
1627 # expressions involving only numbers, the log function, and log(x).
1628 # The series expansion of b**e is computed as follows:
1629 # 1) We express b as f*(1 + g) where f is the leading term of b.
1630 # g has order O(x**d) where d is strictly positive.
1631 # 2) Then b**e = (f**e)*((1 + g)**e).
1632 # (1 + g)**e is computed using binomial series.
1633 from sympy.functions.elementary.exponential import exp, log
1634 from sympy.series.limits import limit
1635 from sympy.series.order import Order
1636 from sympy.core.sympify import sympify
1637 if self.base is S.Exp1:
1638 e_series = self.exp.nseries(x, n=n, logx=logx)
1639 if e_series.is_Order:
1640 return 1 + e_series
1641 e0 = limit(e_series.removeO(), x, 0)
1642 if e0 is S.NegativeInfinity:
1643 return Order(x**n, x)
1644 if e0 is S.Infinity:
1645 return self
1646 t = e_series - e0
1647 exp_series = term = exp(e0)
1648 # series of exp(e0 + t) in t
1649 for i in range(1, n):
1650 term *= t/i
1651 term = term.nseries(x, n=n, logx=logx)
1652 exp_series += term
1653 exp_series += Order(t**n, x)
1654 from sympy.simplify.powsimp import powsimp
1655 return powsimp(exp_series, deep=True, combine='exp')
1656 from sympy.simplify.powsimp import powdenest
1657 from .numbers import _illegal
1658 self = powdenest(self, force=True).trigsimp()
1659 b, e = self.as_base_exp()
1661 if e.has(*_illegal):
1662 raise PoleError()
1664 if e.has(x):
1665 return exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
1667 if logx is not None and b.has(log):
1668 from .symbol import Wild
1669 c, ex = symbols('c, ex', cls=Wild, exclude=[x])
1670 b = b.replace(log(c*x**ex), log(c) + ex*logx)
1671 self = b**e
1673 b = b.removeO()
1674 try:
1675 from sympy.functions.special.gamma_functions import polygamma
1676 if b.has(polygamma, S.EulerGamma) and logx is not None:
1677 raise ValueError()
1678 _, m = b.leadterm(x)
1679 except (ValueError, NotImplementedError, PoleError):
1680 b = b._eval_nseries(x, n=max(2, n), logx=logx, cdir=cdir).removeO()
1681 if b.has(S.NaN, S.ComplexInfinity):
1682 raise NotImplementedError()
1683 _, m = b.leadterm(x)
1685 if e.has(log):
1686 from sympy.simplify.simplify import logcombine
1687 e = logcombine(e).cancel()
1689 if not (m.is_zero or e.is_number and e.is_real):
1690 if self == self._eval_as_leading_term(x, logx=logx, cdir=cdir):
1691 res = exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
1692 if res == exp(e*log(b)):
1693 return self
1694 return res
1696 f = b.as_leading_term(x, logx=logx)
1697 g = (b/f - S.One).cancel(expand=False)
1698 if not m.is_number:
1699 raise NotImplementedError()
1700 maxpow = n - m*e
1701 if maxpow.has(Symbol):
1702 maxpow = sympify(n)
1704 if maxpow.is_negative:
1705 return Order(x**(m*e), x)
1707 if g.is_zero:
1708 r = f**e
1709 if r != self:
1710 r += Order(x**n, x)
1711 return r
1713 def coeff_exp(term, x):
1714 coeff, exp = S.One, S.Zero
1715 for factor in Mul.make_args(term):
1716 if factor.has(x):
1717 base, exp = factor.as_base_exp()
1718 if base != x:
1719 try:
1720 return term.leadterm(x)
1721 except ValueError:
1722 return term, S.Zero
1723 else:
1724 coeff *= factor
1725 return coeff, exp
1727 def mul(d1, d2):
1728 res = {}
1729 for e1, e2 in product(d1, d2):
1730 ex = e1 + e2
1731 if ex < maxpow:
1732 res[ex] = res.get(ex, S.Zero) + d1[e1]*d2[e2]
1733 return res
1735 try:
1736 c, d = g.leadterm(x, logx=logx)
1737 except (ValueError, NotImplementedError):
1738 if limit(g/x**maxpow, x, 0) == 0:
1739 # g has higher order zero
1740 return f**e + e*f**e*g # first term of binomial series
1741 else:
1742 raise NotImplementedError()
1743 if c.is_Float and d == S.Zero:
1744 # Convert floats like 0.5 to exact SymPy numbers like S.Half, to
1745 # prevent rounding errors which can induce wrong values of d leading
1746 # to a NotImplementedError being returned from the block below.
1747 from sympy.simplify.simplify import nsimplify
1748 _, d = nsimplify(g).leadterm(x, logx=logx)
1749 if not d.is_positive:
1750 g = g.simplify()
1751 if g.is_zero:
1752 return f**e
1753 _, d = g.leadterm(x, logx=logx)
1754 if not d.is_positive:
1755 g = ((b - f)/f).expand()
1756 _, d = g.leadterm(x, logx=logx)
1757 if not d.is_positive:
1758 raise NotImplementedError()
1760 from sympy.functions.elementary.integers import ceiling
1761 gpoly = g._eval_nseries(x, n=ceiling(maxpow), logx=logx, cdir=cdir).removeO()
1762 gterms = {}
1764 for term in Add.make_args(gpoly):
1765 co1, e1 = coeff_exp(term, x)
1766 gterms[e1] = gterms.get(e1, S.Zero) + co1
1768 k = S.One
1769 terms = {S.Zero: S.One}
1770 tk = gterms
1772 from sympy.functions.combinatorial.factorials import factorial, ff
1774 while (k*d - maxpow).is_negative:
1775 coeff = ff(e, k)/factorial(k)
1776 for ex in tk:
1777 terms[ex] = terms.get(ex, S.Zero) + coeff*tk[ex]
1778 tk = mul(tk, gterms)
1779 k += S.One
1781 from sympy.functions.elementary.complexes import im
1783 if not e.is_integer and m.is_zero and f.is_negative:
1784 ndir = (b - f).dir(x, cdir)
1785 if im(ndir).is_negative:
1786 inco, inex = coeff_exp(f**e*(-1)**(-2*e), x)
1787 elif im(ndir).is_zero:
1788 inco, inex = coeff_exp(exp(e*log(b)).as_leading_term(x, logx=logx, cdir=cdir), x)
1789 else:
1790 inco, inex = coeff_exp(f**e, x)
1791 else:
1792 inco, inex = coeff_exp(f**e, x)
1793 res = S.Zero
1795 for e1 in terms:
1796 ex = e1 + inex
1797 res += terms[e1]*inco*x**(ex)
1799 if not (e.is_integer and e.is_positive and (e*d - n).is_nonpositive and
1800 res == _mexpand(self)):
1801 try:
1802 res += Order(x**n, x)
1803 except NotImplementedError:
1804 return exp(e*log(b))._eval_nseries(x, n=n, logx=logx, cdir=cdir)
1805 return res
1807 def _eval_as_leading_term(self, x, logx=None, cdir=0):
1808 from sympy.functions.elementary.exponential import exp, log
1809 e = self.exp
1810 b = self.base
1811 if self.base is S.Exp1:
1812 arg = e.as_leading_term(x, logx=logx)
1813 arg0 = arg.subs(x, 0)
1814 if arg0 is S.NaN:
1815 arg0 = arg.limit(x, 0)
1816 if arg0.is_infinite is False:
1817 return S.Exp1**arg0
1818 raise PoleError("Cannot expand %s around 0" % (self))
1819 elif e.has(x):
1820 lt = exp(e * log(b))
1821 return lt.as_leading_term(x, logx=logx, cdir=cdir)
1822 else:
1823 from sympy.functions.elementary.complexes import im
1824 try:
1825 f = b.as_leading_term(x, logx=logx, cdir=cdir)
1826 except PoleError:
1827 return self
1828 if not e.is_integer and f.is_negative and not f.has(x):
1829 ndir = (b - f).dir(x, cdir)
1830 if im(ndir).is_negative:
1831 # Normally, f**e would evaluate to exp(e*log(f)) but on branch cuts
1832 # an other value is expected through the following computation
1833 # exp(e*(log(f) - 2*pi*I)) == f**e*exp(-2*e*pi*I) == f**e*(-1)**(-2*e).
1834 return self.func(f, e) * (-1)**(-2*e)
1835 elif im(ndir).is_zero:
1836 log_leadterm = log(b)._eval_as_leading_term(x, logx=logx, cdir=cdir)
1837 if log_leadterm.is_infinite is False:
1838 return exp(e*log_leadterm)
1839 return self.func(f, e)
1841 @cacheit
1842 def _taylor_term(self, n, x, *previous_terms): # of (1 + x)**e
1843 from sympy.functions.combinatorial.factorials import binomial
1844 return binomial(self.exp, n) * self.func(x, n)
1846 def taylor_term(self, n, x, *previous_terms):
1847 if self.base is not S.Exp1:
1848 return super().taylor_term(n, x, *previous_terms)
1849 if n < 0:
1850 return S.Zero
1851 if n == 0:
1852 return S.One
1853 from .sympify import sympify
1854 x = sympify(x)
1855 if previous_terms:
1856 p = previous_terms[-1]
1857 if p is not None:
1858 return p * x / n
1859 from sympy.functions.combinatorial.factorials import factorial
1860 return x**n/factorial(n)
1862 def _eval_rewrite_as_sin(self, base, exp):
1863 if self.base is S.Exp1:
1864 from sympy.functions.elementary.trigonometric import sin
1865 return sin(S.ImaginaryUnit*self.exp + S.Pi/2) - S.ImaginaryUnit*sin(S.ImaginaryUnit*self.exp)
1867 def _eval_rewrite_as_cos(self, base, exp):
1868 if self.base is S.Exp1:
1869 from sympy.functions.elementary.trigonometric import cos
1870 return cos(S.ImaginaryUnit*self.exp) + S.ImaginaryUnit*cos(S.ImaginaryUnit*self.exp + S.Pi/2)
1872 def _eval_rewrite_as_tanh(self, base, exp):
1873 if self.base is S.Exp1:
1874 from sympy.functions.elementary.hyperbolic import tanh
1875 return (1 + tanh(self.exp/2))/(1 - tanh(self.exp/2))
1877 def _eval_rewrite_as_sqrt(self, base, exp, **kwargs):
1878 from sympy.functions.elementary.trigonometric import sin, cos
1879 if base is not S.Exp1:
1880 return None
1881 if exp.is_Mul:
1882 coeff = exp.coeff(S.Pi * S.ImaginaryUnit)
1883 if coeff and coeff.is_number:
1884 cosine, sine = cos(S.Pi*coeff), sin(S.Pi*coeff)
1885 if not isinstance(cosine, cos) and not isinstance (sine, sin):
1886 return cosine + S.ImaginaryUnit*sine
1888 def as_content_primitive(self, radical=False, clear=True):
1889 """Return the tuple (R, self/R) where R is the positive Rational
1890 extracted from self.
1892 Examples
1893 ========
1895 >>> from sympy import sqrt
1896 >>> sqrt(4 + 4*sqrt(2)).as_content_primitive()
1897 (2, sqrt(1 + sqrt(2)))
1898 >>> sqrt(3 + 3*sqrt(2)).as_content_primitive()
1899 (1, sqrt(3)*sqrt(1 + sqrt(2)))
1901 >>> from sympy import expand_power_base, powsimp, Mul
1902 >>> from sympy.abc import x, y
1904 >>> ((2*x + 2)**2).as_content_primitive()
1905 (4, (x + 1)**2)
1906 >>> (4**((1 + y)/2)).as_content_primitive()
1907 (2, 4**(y/2))
1908 >>> (3**((1 + y)/2)).as_content_primitive()
1909 (1, 3**((y + 1)/2))
1910 >>> (3**((5 + y)/2)).as_content_primitive()
1911 (9, 3**((y + 1)/2))
1912 >>> eq = 3**(2 + 2*x)
1913 >>> powsimp(eq) == eq
1914 True
1915 >>> eq.as_content_primitive()
1916 (9, 3**(2*x))
1917 >>> powsimp(Mul(*_))
1918 3**(2*x + 2)
1920 >>> eq = (2 + 2*x)**y
1921 >>> s = expand_power_base(eq); s.is_Mul, s
1922 (False, (2*x + 2)**y)
1923 >>> eq.as_content_primitive()
1924 (1, (2*(x + 1))**y)
1925 >>> s = expand_power_base(_[1]); s.is_Mul, s
1926 (True, 2**y*(x + 1)**y)
1928 See docstring of Expr.as_content_primitive for more examples.
1929 """
1931 b, e = self.as_base_exp()
1932 b = _keep_coeff(*b.as_content_primitive(radical=radical, clear=clear))
1933 ce, pe = e.as_content_primitive(radical=radical, clear=clear)
1934 if b.is_Rational:
1935 #e
1936 #= ce*pe
1937 #= ce*(h + t)
1938 #= ce*h + ce*t
1939 #=> self
1940 #= b**(ce*h)*b**(ce*t)
1941 #= b**(cehp/cehq)*b**(ce*t)
1942 #= b**(iceh + r/cehq)*b**(ce*t)
1943 #= b**(iceh)*b**(r/cehq)*b**(ce*t)
1944 #= b**(iceh)*b**(ce*t + r/cehq)
1945 h, t = pe.as_coeff_Add()
1946 if h.is_Rational and b != S.Zero:
1947 ceh = ce*h
1948 c = self.func(b, ceh)
1949 r = S.Zero
1950 if not c.is_Rational:
1951 iceh, r = divmod(ceh.p, ceh.q)
1952 c = self.func(b, iceh)
1953 return c, self.func(b, _keep_coeff(ce, t + r/ce/ceh.q))
1954 e = _keep_coeff(ce, pe)
1955 # b**e = (h*t)**e = h**e*t**e = c*m*t**e
1956 if e.is_Rational and b.is_Mul:
1957 h, t = b.as_content_primitive(radical=radical, clear=clear) # h is positive
1958 c, m = self.func(h, e).as_coeff_Mul() # so c is positive
1959 m, me = m.as_base_exp()
1960 if m is S.One or me == e: # probably always true
1961 # return the following, not return c, m*Pow(t, e)
1962 # which would change Pow into Mul; we let SymPy
1963 # decide what to do by using the unevaluated Mul, e.g
1964 # should it stay as sqrt(2 + 2*sqrt(5)) or become
1965 # sqrt(2)*sqrt(1 + sqrt(5))
1966 return c, self.func(_keep_coeff(m, t), e)
1967 return S.One, self.func(b, e)
1969 def is_constant(self, *wrt, **flags):
1970 expr = self
1971 if flags.get('simplify', True):
1972 expr = expr.simplify()
1973 b, e = expr.as_base_exp()
1974 bz = b.equals(0)
1975 if bz: # recalculate with assumptions in case it's unevaluated
1976 new = b**e
1977 if new != expr:
1978 return new.is_constant()
1979 econ = e.is_constant(*wrt)
1980 bcon = b.is_constant(*wrt)
1981 if bcon:
1982 if econ:
1983 return True
1984 bz = b.equals(0)
1985 if bz is False:
1986 return False
1987 elif bcon is None:
1988 return None
1990 return e.equals(0)
1992 def _eval_difference_delta(self, n, step):
1993 b, e = self.args
1994 if e.has(n) and not b.has(n):
1995 new_e = e.subs(n, n + step)
1996 return (b**(new_e - e) - 1) * self
1998power = Dispatcher('power')
1999power.add((object, object), Pow)
2001from .add import Add
2002from .numbers import Integer
2003from .mul import Mul, _keep_coeff
2004from .symbol import Symbol, Dummy, symbols