Coverage for /usr/lib/python3/dist-packages/sympy/core/exprtools.py: 37%
779 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Tools for manipulating of large commutative expressions. """
3from .add import Add
4from .mul import Mul, _keep_coeff
5from .power import Pow
6from .basic import Basic
7from .expr import Expr
8from .function import expand_power_exp
9from .sympify import sympify
10from .numbers import Rational, Integer, Number, I, equal_valued
11from .singleton import S
12from .sorting import default_sort_key, ordered
13from .symbol import Dummy
14from .traversal import preorder_traversal
15from .coreerrors import NonCommutativeExpression
16from .containers import Tuple, Dict
17from sympy.external.gmpy import SYMPY_INTS
18from sympy.utilities.iterables import (common_prefix, common_suffix,
19 variations, iterable, is_sequence)
21from collections import defaultdict
22from typing import Tuple as tTuple
25_eps = Dummy(positive=True)
28def _isnumber(i):
29 return isinstance(i, (SYMPY_INTS, float)) or i.is_Number
32def _monotonic_sign(self):
33 """Return the value closest to 0 that ``self`` may have if all symbols
34 are signed and the result is uniformly the same sign for all values of symbols.
35 If a symbol is only signed but not known to be an
36 integer or the result is 0 then a symbol representative of the sign of self
37 will be returned. Otherwise, None is returned if a) the sign could be positive
38 or negative or b) self is not in one of the following forms:
40 - L(x, y, ...) + A: a function linear in all symbols x, y, ... with an
41 additive constant; if A is zero then the function can be a monomial whose
42 sign is monotonic over the range of the variables, e.g. (x + 1)**3 if x is
43 nonnegative.
44 - A/L(x, y, ...) + B: the inverse of a function linear in all symbols x, y, ...
45 that does not have a sign change from positive to negative for any set
46 of values for the variables.
47 - M(x, y, ...) + A: a monomial M whose factors are all signed and a constant, A.
48 - A/M(x, y, ...) + B: the inverse of a monomial and constants A and B.
49 - P(x): a univariate polynomial
51 Examples
52 ========
54 >>> from sympy.core.exprtools import _monotonic_sign as F
55 >>> from sympy import Dummy
56 >>> nn = Dummy(integer=True, nonnegative=True)
57 >>> p = Dummy(integer=True, positive=True)
58 >>> p2 = Dummy(integer=True, positive=True)
59 >>> F(nn + 1)
60 1
61 >>> F(p - 1)
62 _nneg
63 >>> F(nn*p + 1)
64 1
65 >>> F(p2*p + 1)
66 2
67 >>> F(nn - 1) # could be negative, zero or positive
68 """
69 if not self.is_extended_real:
70 return
72 if (-self).is_Symbol:
73 rv = _monotonic_sign(-self)
74 return rv if rv is None else -rv
76 if not self.is_Add and self.as_numer_denom()[1].is_number:
77 s = self
78 if s.is_prime:
79 if s.is_odd:
80 return Integer(3)
81 else:
82 return Integer(2)
83 elif s.is_composite:
84 if s.is_odd:
85 return Integer(9)
86 else:
87 return Integer(4)
88 elif s.is_positive:
89 if s.is_even:
90 if s.is_prime is False:
91 return Integer(4)
92 else:
93 return Integer(2)
94 elif s.is_integer:
95 return S.One
96 else:
97 return _eps
98 elif s.is_extended_negative:
99 if s.is_even:
100 return Integer(-2)
101 elif s.is_integer:
102 return S.NegativeOne
103 else:
104 return -_eps
105 if s.is_zero or s.is_extended_nonpositive or s.is_extended_nonnegative:
106 return S.Zero
107 return None
109 # univariate polynomial
110 free = self.free_symbols
111 if len(free) == 1:
112 if self.is_polynomial():
113 from sympy.polys.polytools import real_roots
114 from sympy.polys.polyroots import roots
115 from sympy.polys.polyerrors import PolynomialError
116 x = free.pop()
117 x0 = _monotonic_sign(x)
118 if x0 in (_eps, -_eps):
119 x0 = S.Zero
120 if x0 is not None:
121 d = self.diff(x)
122 if d.is_number:
123 currentroots = []
124 else:
125 try:
126 currentroots = real_roots(d)
127 except (PolynomialError, NotImplementedError):
128 currentroots = [r for r in roots(d, x) if r.is_extended_real]
129 y = self.subs(x, x0)
130 if x.is_nonnegative and all(
131 (r - x0).is_nonpositive for r in currentroots):
132 if y.is_nonnegative and d.is_positive:
133 if y:
134 return y if y.is_positive else Dummy('pos', positive=True)
135 else:
136 return Dummy('nneg', nonnegative=True)
137 if y.is_nonpositive and d.is_negative:
138 if y:
139 return y if y.is_negative else Dummy('neg', negative=True)
140 else:
141 return Dummy('npos', nonpositive=True)
142 elif x.is_nonpositive and all(
143 (r - x0).is_nonnegative for r in currentroots):
144 if y.is_nonnegative and d.is_negative:
145 if y:
146 return Dummy('pos', positive=True)
147 else:
148 return Dummy('nneg', nonnegative=True)
149 if y.is_nonpositive and d.is_positive:
150 if y:
151 return Dummy('neg', negative=True)
152 else:
153 return Dummy('npos', nonpositive=True)
154 else:
155 n, d = self.as_numer_denom()
156 den = None
157 if n.is_number:
158 den = _monotonic_sign(d)
159 elif not d.is_number:
160 if _monotonic_sign(n) is not None:
161 den = _monotonic_sign(d)
162 if den is not None and (den.is_positive or den.is_negative):
163 v = n*den
164 if v.is_positive:
165 return Dummy('pos', positive=True)
166 elif v.is_nonnegative:
167 return Dummy('nneg', nonnegative=True)
168 elif v.is_negative:
169 return Dummy('neg', negative=True)
170 elif v.is_nonpositive:
171 return Dummy('npos', nonpositive=True)
172 return None
174 # multivariate
175 c, a = self.as_coeff_Add()
176 v = None
177 if not a.is_polynomial():
178 # F/A or A/F where A is a number and F is a signed, rational monomial
179 n, d = a.as_numer_denom()
180 if not (n.is_number or d.is_number):
181 return
182 if (
183 a.is_Mul or a.is_Pow) and \
184 a.is_rational and \
185 all(p.exp.is_Integer for p in a.atoms(Pow) if p.is_Pow) and \
186 (a.is_positive or a.is_negative):
187 v = S.One
188 for ai in Mul.make_args(a):
189 if ai.is_number:
190 v *= ai
191 continue
192 reps = {}
193 for x in ai.free_symbols:
194 reps[x] = _monotonic_sign(x)
195 if reps[x] is None:
196 return
197 v *= ai.subs(reps)
198 elif c:
199 # signed linear expression
200 if not any(p for p in a.atoms(Pow) if not p.is_number) and (a.is_nonpositive or a.is_nonnegative):
201 free = list(a.free_symbols)
202 p = {}
203 for i in free:
204 v = _monotonic_sign(i)
205 if v is None:
206 return
207 p[i] = v or (_eps if i.is_nonnegative else -_eps)
208 v = a.xreplace(p)
209 if v is not None:
210 rv = v + c
211 if v.is_nonnegative and rv.is_positive:
212 return rv.subs(_eps, 0)
213 if v.is_nonpositive and rv.is_negative:
214 return rv.subs(_eps, 0)
217def decompose_power(expr: Expr) -> tTuple[Expr, int]:
218 """
219 Decompose power into symbolic base and integer exponent.
221 Examples
222 ========
224 >>> from sympy.core.exprtools import decompose_power
225 >>> from sympy.abc import x, y
226 >>> from sympy import exp
228 >>> decompose_power(x)
229 (x, 1)
230 >>> decompose_power(x**2)
231 (x, 2)
232 >>> decompose_power(exp(2*y/3))
233 (exp(y/3), 2)
235 """
236 base, exp = expr.as_base_exp()
238 if exp.is_Number:
239 if exp.is_Rational:
240 if not exp.is_Integer:
241 base = Pow(base, Rational(1, exp.q)) # type: ignore
242 e = exp.p # type: ignore
243 else:
244 base, e = expr, 1
245 else:
246 exp, tail = exp.as_coeff_Mul(rational=True)
248 if exp is S.NegativeOne:
249 base, e = Pow(base, tail), -1
250 elif exp is not S.One:
251 # todo: after dropping python 3.7 support, use overload and Literal
252 # in as_coeff_Mul to make exp Rational, and remove these 2 ignores
253 tail = _keep_coeff(Rational(1, exp.q), tail) # type: ignore
254 base, e = Pow(base, tail), exp.p # type: ignore
255 else:
256 base, e = expr, 1
258 return base, e
261def decompose_power_rat(expr: Expr) -> tTuple[Expr, Rational]:
262 """
263 Decompose power into symbolic base and rational exponent;
264 if the exponent is not a Rational, then separate only the
265 integer coefficient.
267 Examples
268 ========
270 >>> from sympy.core.exprtools import decompose_power_rat
271 >>> from sympy.abc import x
272 >>> from sympy import sqrt, exp
274 >>> decompose_power_rat(sqrt(x))
275 (x, 1/2)
276 >>> decompose_power_rat(exp(-3*x/2))
277 (exp(x/2), -3)
279 """
280 _ = base, exp = expr.as_base_exp()
281 return _ if exp.is_Rational else decompose_power(expr)
284class Factors:
285 """Efficient representation of ``f_1*f_2*...*f_n``."""
287 __slots__ = ('factors', 'gens')
289 def __init__(self, factors=None): # Factors
290 """Initialize Factors from dict or expr.
292 Examples
293 ========
295 >>> from sympy.core.exprtools import Factors
296 >>> from sympy.abc import x
297 >>> from sympy import I
298 >>> e = 2*x**3
299 >>> Factors(e)
300 Factors({2: 1, x: 3})
301 >>> Factors(e.as_powers_dict())
302 Factors({2: 1, x: 3})
303 >>> f = _
304 >>> f.factors # underlying dictionary
305 {2: 1, x: 3}
306 >>> f.gens # base of each factor
307 frozenset({2, x})
308 >>> Factors(0)
309 Factors({0: 1})
310 >>> Factors(I)
311 Factors({I: 1})
313 Notes
314 =====
316 Although a dictionary can be passed, only minimal checking is
317 performed: powers of -1 and I are made canonical.
319 """
320 if isinstance(factors, (SYMPY_INTS, float)):
321 factors = S(factors)
322 if isinstance(factors, Factors):
323 factors = factors.factors.copy()
324 elif factors in (None, S.One):
325 factors = {}
326 elif factors is S.Zero or factors == 0:
327 factors = {S.Zero: S.One}
328 elif isinstance(factors, Number):
329 n = factors
330 factors = {}
331 if n < 0:
332 factors[S.NegativeOne] = S.One
333 n = -n
334 if n is not S.One:
335 if n.is_Float or n.is_Integer or n is S.Infinity:
336 factors[n] = S.One
337 elif n.is_Rational:
338 # since we're processing Numbers, the denominator is
339 # stored with a negative exponent; all other factors
340 # are left .
341 if n.p != 1:
342 factors[Integer(n.p)] = S.One
343 factors[Integer(n.q)] = S.NegativeOne
344 else:
345 raise ValueError('Expected Float|Rational|Integer, not %s' % n)
346 elif isinstance(factors, Basic) and not factors.args:
347 factors = {factors: S.One}
348 elif isinstance(factors, Expr):
349 c, nc = factors.args_cnc()
350 i = c.count(I)
351 for _ in range(i):
352 c.remove(I)
353 factors = dict(Mul._from_args(c).as_powers_dict())
354 # Handle all rational Coefficients
355 for f in list(factors.keys()):
356 if isinstance(f, Rational) and not isinstance(f, Integer):
357 p, q = Integer(f.p), Integer(f.q)
358 factors[p] = (factors[p] if p in factors else S.Zero) + factors[f]
359 factors[q] = (factors[q] if q in factors else S.Zero) - factors[f]
360 factors.pop(f)
361 if i:
362 factors[I] = factors.get(I, S.Zero) + i
363 if nc:
364 factors[Mul(*nc, evaluate=False)] = S.One
365 else:
366 factors = factors.copy() # /!\ should be dict-like
368 # tidy up -/+1 and I exponents if Rational
370 handle = [k for k in factors if k is I or k in (-1, 1)]
371 if handle:
372 i1 = S.One
373 for k in handle:
374 if not _isnumber(factors[k]):
375 continue
376 i1 *= k**factors.pop(k)
377 if i1 is not S.One:
378 for a in i1.args if i1.is_Mul else [i1]: # at worst, -1.0*I*(-1)**e
379 if a is S.NegativeOne:
380 factors[a] = S.One
381 elif a is I:
382 factors[I] = S.One
383 elif a.is_Pow:
384 factors[a.base] = factors.get(a.base, S.Zero) + a.exp
385 elif equal_valued(a, 1):
386 factors[a] = S.One
387 elif equal_valued(a, -1):
388 factors[-a] = S.One
389 factors[S.NegativeOne] = S.One
390 else:
391 raise ValueError('unexpected factor in i1: %s' % a)
393 self.factors = factors
394 keys = getattr(factors, 'keys', None)
395 if keys is None:
396 raise TypeError('expecting Expr or dictionary')
397 self.gens = frozenset(keys())
399 def __hash__(self): # Factors
400 keys = tuple(ordered(self.factors.keys()))
401 values = [self.factors[k] for k in keys]
402 return hash((keys, values))
404 def __repr__(self): # Factors
405 return "Factors({%s})" % ', '.join(
406 ['%s: %s' % (k, v) for k, v in ordered(self.factors.items())])
408 @property
409 def is_zero(self): # Factors
410 """
411 >>> from sympy.core.exprtools import Factors
412 >>> Factors(0).is_zero
413 True
414 """
415 f = self.factors
416 return len(f) == 1 and S.Zero in f
418 @property
419 def is_one(self): # Factors
420 """
421 >>> from sympy.core.exprtools import Factors
422 >>> Factors(1).is_one
423 True
424 """
425 return not self.factors
427 def as_expr(self): # Factors
428 """Return the underlying expression.
430 Examples
431 ========
433 >>> from sympy.core.exprtools import Factors
434 >>> from sympy.abc import x, y
435 >>> Factors((x*y**2).as_powers_dict()).as_expr()
436 x*y**2
438 """
440 args = []
441 for factor, exp in self.factors.items():
442 if exp != 1:
443 if isinstance(exp, Integer):
444 b, e = factor.as_base_exp()
445 e = _keep_coeff(exp, e)
446 args.append(b**e)
447 else:
448 args.append(factor**exp)
449 else:
450 args.append(factor)
451 return Mul(*args)
453 def mul(self, other): # Factors
454 """Return Factors of ``self * other``.
456 Examples
457 ========
459 >>> from sympy.core.exprtools import Factors
460 >>> from sympy.abc import x, y, z
461 >>> a = Factors((x*y**2).as_powers_dict())
462 >>> b = Factors((x*y/z).as_powers_dict())
463 >>> a.mul(b)
464 Factors({x: 2, y: 3, z: -1})
465 >>> a*b
466 Factors({x: 2, y: 3, z: -1})
467 """
468 if not isinstance(other, Factors):
469 other = Factors(other)
470 if any(f.is_zero for f in (self, other)):
471 return Factors(S.Zero)
472 factors = dict(self.factors)
474 for factor, exp in other.factors.items():
475 if factor in factors:
476 exp = factors[factor] + exp
478 if not exp:
479 del factors[factor]
480 continue
482 factors[factor] = exp
484 return Factors(factors)
486 def normal(self, other):
487 """Return ``self`` and ``other`` with ``gcd`` removed from each.
488 The only differences between this and method ``div`` is that this
489 is 1) optimized for the case when there are few factors in common and
490 2) this does not raise an error if ``other`` is zero.
492 See Also
493 ========
494 div
496 """
497 if not isinstance(other, Factors):
498 other = Factors(other)
499 if other.is_zero:
500 return (Factors(), Factors(S.Zero))
501 if self.is_zero:
502 return (Factors(S.Zero), Factors())
504 self_factors = dict(self.factors)
505 other_factors = dict(other.factors)
507 for factor, self_exp in self.factors.items():
508 try:
509 other_exp = other.factors[factor]
510 except KeyError:
511 continue
513 exp = self_exp - other_exp
515 if not exp:
516 del self_factors[factor]
517 del other_factors[factor]
518 elif _isnumber(exp):
519 if exp > 0:
520 self_factors[factor] = exp
521 del other_factors[factor]
522 else:
523 del self_factors[factor]
524 other_factors[factor] = -exp
525 else:
526 r = self_exp.extract_additively(other_exp)
527 if r is not None:
528 if r:
529 self_factors[factor] = r
530 del other_factors[factor]
531 else: # should be handled already
532 del self_factors[factor]
533 del other_factors[factor]
534 else:
535 sc, sa = self_exp.as_coeff_Add()
536 if sc:
537 oc, oa = other_exp.as_coeff_Add()
538 diff = sc - oc
539 if diff > 0:
540 self_factors[factor] -= oc
541 other_exp = oa
542 elif diff < 0:
543 self_factors[factor] -= sc
544 other_factors[factor] -= sc
545 other_exp = oa - diff
546 else:
547 self_factors[factor] = sa
548 other_exp = oa
549 if other_exp:
550 other_factors[factor] = other_exp
551 else:
552 del other_factors[factor]
554 return Factors(self_factors), Factors(other_factors)
556 def div(self, other): # Factors
557 """Return ``self`` and ``other`` with ``gcd`` removed from each.
558 This is optimized for the case when there are many factors in common.
560 Examples
561 ========
563 >>> from sympy.core.exprtools import Factors
564 >>> from sympy.abc import x, y, z
565 >>> from sympy import S
567 >>> a = Factors((x*y**2).as_powers_dict())
568 >>> a.div(a)
569 (Factors({}), Factors({}))
570 >>> a.div(x*z)
571 (Factors({y: 2}), Factors({z: 1}))
573 The ``/`` operator only gives ``quo``:
575 >>> a/x
576 Factors({y: 2})
578 Factors treats its factors as though they are all in the numerator, so
579 if you violate this assumption the results will be correct but will
580 not strictly correspond to the numerator and denominator of the ratio:
582 >>> a.div(x/z)
583 (Factors({y: 2}), Factors({z: -1}))
585 Factors is also naive about bases: it does not attempt any denesting
586 of Rational-base terms, for example the following does not become
587 2**(2*x)/2.
589 >>> Factors(2**(2*x + 2)).div(S(8))
590 (Factors({2: 2*x + 2}), Factors({8: 1}))
592 factor_terms can clean up such Rational-bases powers:
594 >>> from sympy import factor_terms
595 >>> n, d = Factors(2**(2*x + 2)).div(S(8))
596 >>> n.as_expr()/d.as_expr()
597 2**(2*x + 2)/8
598 >>> factor_terms(_)
599 2**(2*x)/2
601 """
602 quo, rem = dict(self.factors), {}
604 if not isinstance(other, Factors):
605 other = Factors(other)
606 if other.is_zero:
607 raise ZeroDivisionError
608 if self.is_zero:
609 return (Factors(S.Zero), Factors())
611 for factor, exp in other.factors.items():
612 if factor in quo:
613 d = quo[factor] - exp
614 if _isnumber(d):
615 if d <= 0:
616 del quo[factor]
618 if d >= 0:
619 if d:
620 quo[factor] = d
622 continue
624 exp = -d
626 else:
627 r = quo[factor].extract_additively(exp)
628 if r is not None:
629 if r:
630 quo[factor] = r
631 else: # should be handled already
632 del quo[factor]
633 else:
634 other_exp = exp
635 sc, sa = quo[factor].as_coeff_Add()
636 if sc:
637 oc, oa = other_exp.as_coeff_Add()
638 diff = sc - oc
639 if diff > 0:
640 quo[factor] -= oc
641 other_exp = oa
642 elif diff < 0:
643 quo[factor] -= sc
644 other_exp = oa - diff
645 else:
646 quo[factor] = sa
647 other_exp = oa
648 if other_exp:
649 rem[factor] = other_exp
650 else:
651 assert factor not in rem
652 continue
654 rem[factor] = exp
656 return Factors(quo), Factors(rem)
658 def quo(self, other): # Factors
659 """Return numerator Factor of ``self / other``.
661 Examples
662 ========
664 >>> from sympy.core.exprtools import Factors
665 >>> from sympy.abc import x, y, z
666 >>> a = Factors((x*y**2).as_powers_dict())
667 >>> b = Factors((x*y/z).as_powers_dict())
668 >>> a.quo(b) # same as a/b
669 Factors({y: 1})
670 """
671 return self.div(other)[0]
673 def rem(self, other): # Factors
674 """Return denominator Factors of ``self / other``.
676 Examples
677 ========
679 >>> from sympy.core.exprtools import Factors
680 >>> from sympy.abc import x, y, z
681 >>> a = Factors((x*y**2).as_powers_dict())
682 >>> b = Factors((x*y/z).as_powers_dict())
683 >>> a.rem(b)
684 Factors({z: -1})
685 >>> a.rem(a)
686 Factors({})
687 """
688 return self.div(other)[1]
690 def pow(self, other): # Factors
691 """Return self raised to a non-negative integer power.
693 Examples
694 ========
696 >>> from sympy.core.exprtools import Factors
697 >>> from sympy.abc import x, y
698 >>> a = Factors((x*y**2).as_powers_dict())
699 >>> a**2
700 Factors({x: 2, y: 4})
702 """
703 if isinstance(other, Factors):
704 other = other.as_expr()
705 if other.is_Integer:
706 other = int(other)
707 if isinstance(other, SYMPY_INTS) and other >= 0:
708 factors = {}
710 if other:
711 for factor, exp in self.factors.items():
712 factors[factor] = exp*other
714 return Factors(factors)
715 else:
716 raise ValueError("expected non-negative integer, got %s" % other)
718 def gcd(self, other): # Factors
719 """Return Factors of ``gcd(self, other)``. The keys are
720 the intersection of factors with the minimum exponent for
721 each factor.
723 Examples
724 ========
726 >>> from sympy.core.exprtools import Factors
727 >>> from sympy.abc import x, y, z
728 >>> a = Factors((x*y**2).as_powers_dict())
729 >>> b = Factors((x*y/z).as_powers_dict())
730 >>> a.gcd(b)
731 Factors({x: 1, y: 1})
732 """
733 if not isinstance(other, Factors):
734 other = Factors(other)
735 if other.is_zero:
736 return Factors(self.factors)
738 factors = {}
740 for factor, exp in self.factors.items():
741 factor, exp = sympify(factor), sympify(exp)
742 if factor in other.factors:
743 lt = (exp - other.factors[factor]).is_negative
744 if lt == True:
745 factors[factor] = exp
746 elif lt == False:
747 factors[factor] = other.factors[factor]
749 return Factors(factors)
751 def lcm(self, other): # Factors
752 """Return Factors of ``lcm(self, other)`` which are
753 the union of factors with the maximum exponent for
754 each factor.
756 Examples
757 ========
759 >>> from sympy.core.exprtools import Factors
760 >>> from sympy.abc import x, y, z
761 >>> a = Factors((x*y**2).as_powers_dict())
762 >>> b = Factors((x*y/z).as_powers_dict())
763 >>> a.lcm(b)
764 Factors({x: 1, y: 2, z: -1})
765 """
766 if not isinstance(other, Factors):
767 other = Factors(other)
768 if any(f.is_zero for f in (self, other)):
769 return Factors(S.Zero)
771 factors = dict(self.factors)
773 for factor, exp in other.factors.items():
774 if factor in factors:
775 exp = max(exp, factors[factor])
777 factors[factor] = exp
779 return Factors(factors)
781 def __mul__(self, other): # Factors
782 return self.mul(other)
784 def __divmod__(self, other): # Factors
785 return self.div(other)
787 def __truediv__(self, other): # Factors
788 return self.quo(other)
790 def __mod__(self, other): # Factors
791 return self.rem(other)
793 def __pow__(self, other): # Factors
794 return self.pow(other)
796 def __eq__(self, other): # Factors
797 if not isinstance(other, Factors):
798 other = Factors(other)
799 return self.factors == other.factors
801 def __ne__(self, other): # Factors
802 return not self == other
805class Term:
806 """Efficient representation of ``coeff*(numer/denom)``. """
808 __slots__ = ('coeff', 'numer', 'denom')
810 def __init__(self, term, numer=None, denom=None): # Term
811 if numer is None and denom is None:
812 if not term.is_commutative:
813 raise NonCommutativeExpression(
814 'commutative expression expected')
816 coeff, factors = term.as_coeff_mul()
817 numer, denom = defaultdict(int), defaultdict(int)
819 for factor in factors:
820 base, exp = decompose_power(factor)
822 if base.is_Add:
823 cont, base = base.primitive()
824 coeff *= cont**exp
826 if exp > 0:
827 numer[base] += exp
828 else:
829 denom[base] += -exp
831 numer = Factors(numer)
832 denom = Factors(denom)
833 else:
834 coeff = term
836 if numer is None:
837 numer = Factors()
839 if denom is None:
840 denom = Factors()
842 self.coeff = coeff
843 self.numer = numer
844 self.denom = denom
846 def __hash__(self): # Term
847 return hash((self.coeff, self.numer, self.denom))
849 def __repr__(self): # Term
850 return "Term(%s, %s, %s)" % (self.coeff, self.numer, self.denom)
852 def as_expr(self): # Term
853 return self.coeff*(self.numer.as_expr()/self.denom.as_expr())
855 def mul(self, other): # Term
856 coeff = self.coeff*other.coeff
857 numer = self.numer.mul(other.numer)
858 denom = self.denom.mul(other.denom)
860 numer, denom = numer.normal(denom)
862 return Term(coeff, numer, denom)
864 def inv(self): # Term
865 return Term(1/self.coeff, self.denom, self.numer)
867 def quo(self, other): # Term
868 return self.mul(other.inv())
870 def pow(self, other): # Term
871 if other < 0:
872 return self.inv().pow(-other)
873 else:
874 return Term(self.coeff ** other,
875 self.numer.pow(other),
876 self.denom.pow(other))
878 def gcd(self, other): # Term
879 return Term(self.coeff.gcd(other.coeff),
880 self.numer.gcd(other.numer),
881 self.denom.gcd(other.denom))
883 def lcm(self, other): # Term
884 return Term(self.coeff.lcm(other.coeff),
885 self.numer.lcm(other.numer),
886 self.denom.lcm(other.denom))
888 def __mul__(self, other): # Term
889 if isinstance(other, Term):
890 return self.mul(other)
891 else:
892 return NotImplemented
894 def __truediv__(self, other): # Term
895 if isinstance(other, Term):
896 return self.quo(other)
897 else:
898 return NotImplemented
900 def __pow__(self, other): # Term
901 if isinstance(other, SYMPY_INTS):
902 return self.pow(other)
903 else:
904 return NotImplemented
906 def __eq__(self, other): # Term
907 return (self.coeff == other.coeff and
908 self.numer == other.numer and
909 self.denom == other.denom)
911 def __ne__(self, other): # Term
912 return not self == other
915def _gcd_terms(terms, isprimitive=False, fraction=True):
916 """Helper function for :func:`gcd_terms`.
918 Parameters
919 ==========
921 isprimitive : boolean, optional
922 If ``isprimitive`` is True then the call to primitive
923 for an Add will be skipped. This is useful when the
924 content has already been extracted.
926 fraction : boolean, optional
927 If ``fraction`` is True then the expression will appear over a common
928 denominator, the lcm of all term denominators.
929 """
931 if isinstance(terms, Basic) and not isinstance(terms, Tuple):
932 terms = Add.make_args(terms)
934 terms = list(map(Term, [t for t in terms if t]))
936 # there is some simplification that may happen if we leave this
937 # here rather than duplicate it before the mapping of Term onto
938 # the terms
939 if len(terms) == 0:
940 return S.Zero, S.Zero, S.One
942 if len(terms) == 1:
943 cont = terms[0].coeff
944 numer = terms[0].numer.as_expr()
945 denom = terms[0].denom.as_expr()
947 else:
948 cont = terms[0]
949 for term in terms[1:]:
950 cont = cont.gcd(term)
952 for i, term in enumerate(terms):
953 terms[i] = term.quo(cont)
955 if fraction:
956 denom = terms[0].denom
958 for term in terms[1:]:
959 denom = denom.lcm(term.denom)
961 numers = []
962 for term in terms:
963 numer = term.numer.mul(denom.quo(term.denom))
964 numers.append(term.coeff*numer.as_expr())
965 else:
966 numers = [t.as_expr() for t in terms]
967 denom = Term(S.One).numer
969 cont = cont.as_expr()
970 numer = Add(*numers)
971 denom = denom.as_expr()
973 if not isprimitive and numer.is_Add:
974 _cont, numer = numer.primitive()
975 cont *= _cont
977 return cont, numer, denom
980def gcd_terms(terms, isprimitive=False, clear=True, fraction=True):
981 """Compute the GCD of ``terms`` and put them together.
983 Parameters
984 ==========
986 terms : Expr
987 Can be an expression or a non-Basic sequence of expressions
988 which will be handled as though they are terms from a sum.
990 isprimitive : bool, optional
991 If ``isprimitive`` is True the _gcd_terms will not run the primitive
992 method on the terms.
994 clear : bool, optional
995 It controls the removal of integers from the denominator of an Add
996 expression. When True (default), all numerical denominator will be cleared;
997 when False the denominators will be cleared only if all terms had numerical
998 denominators other than 1.
1000 fraction : bool, optional
1001 When True (default), will put the expression over a common
1002 denominator.
1004 Examples
1005 ========
1007 >>> from sympy import gcd_terms
1008 >>> from sympy.abc import x, y
1010 >>> gcd_terms((x + 1)**2*y + (x + 1)*y**2)
1011 y*(x + 1)*(x + y + 1)
1012 >>> gcd_terms(x/2 + 1)
1013 (x + 2)/2
1014 >>> gcd_terms(x/2 + 1, clear=False)
1015 x/2 + 1
1016 >>> gcd_terms(x/2 + y/2, clear=False)
1017 (x + y)/2
1018 >>> gcd_terms(x/2 + 1/x)
1019 (x**2 + 2)/(2*x)
1020 >>> gcd_terms(x/2 + 1/x, fraction=False)
1021 (x + 2/x)/2
1022 >>> gcd_terms(x/2 + 1/x, fraction=False, clear=False)
1023 x/2 + 1/x
1025 >>> gcd_terms(x/2/y + 1/x/y)
1026 (x**2 + 2)/(2*x*y)
1027 >>> gcd_terms(x/2/y + 1/x/y, clear=False)
1028 (x**2/2 + 1)/(x*y)
1029 >>> gcd_terms(x/2/y + 1/x/y, clear=False, fraction=False)
1030 (x/2 + 1/x)/y
1032 The ``clear`` flag was ignored in this case because the returned
1033 expression was a rational expression, not a simple sum.
1035 See Also
1036 ========
1038 factor_terms, sympy.polys.polytools.terms_gcd
1040 """
1041 def mask(terms):
1042 """replace nc portions of each term with a unique Dummy symbols
1043 and return the replacements to restore them"""
1044 args = [(a, []) if a.is_commutative else a.args_cnc() for a in terms]
1045 reps = []
1046 for i, (c, nc) in enumerate(args):
1047 if nc:
1048 nc = Mul(*nc)
1049 d = Dummy()
1050 reps.append((d, nc))
1051 c.append(d)
1052 args[i] = Mul(*c)
1053 else:
1054 args[i] = c
1055 return args, dict(reps)
1057 isadd = isinstance(terms, Add)
1058 addlike = isadd or not isinstance(terms, Basic) and \
1059 is_sequence(terms, include=set) and \
1060 not isinstance(terms, Dict)
1062 if addlike:
1063 if isadd: # i.e. an Add
1064 terms = list(terms.args)
1065 else:
1066 terms = sympify(terms)
1067 terms, reps = mask(terms)
1068 cont, numer, denom = _gcd_terms(terms, isprimitive, fraction)
1069 numer = numer.xreplace(reps)
1070 coeff, factors = cont.as_coeff_Mul()
1071 if not clear:
1072 c, _coeff = coeff.as_coeff_Mul()
1073 if not c.is_Integer and not clear and numer.is_Add:
1074 n, d = c.as_numer_denom()
1075 _numer = numer/d
1076 if any(a.as_coeff_Mul()[0].is_Integer
1077 for a in _numer.args):
1078 numer = _numer
1079 coeff = n*_coeff
1080 return _keep_coeff(coeff, factors*numer/denom, clear=clear)
1082 if not isinstance(terms, Basic):
1083 return terms
1085 if terms.is_Atom:
1086 return terms
1088 if terms.is_Mul:
1089 c, args = terms.as_coeff_mul()
1090 return _keep_coeff(c, Mul(*[gcd_terms(i, isprimitive, clear, fraction)
1091 for i in args]), clear=clear)
1093 def handle(a):
1094 # don't treat internal args like terms of an Add
1095 if not isinstance(a, Expr):
1096 if isinstance(a, Basic):
1097 if not a.args:
1098 return a
1099 return a.func(*[handle(i) for i in a.args])
1100 return type(a)([handle(i) for i in a])
1101 return gcd_terms(a, isprimitive, clear, fraction)
1103 if isinstance(terms, Dict):
1104 return Dict(*[(k, handle(v)) for k, v in terms.args])
1105 return terms.func(*[handle(i) for i in terms.args])
1108def _factor_sum_int(expr, **kwargs):
1109 """Return Sum or Integral object with factors that are not
1110 in the wrt variables removed. In cases where there are additive
1111 terms in the function of the object that are independent, the
1112 object will be separated into two objects.
1114 Examples
1115 ========
1117 >>> from sympy import Sum, factor_terms
1118 >>> from sympy.abc import x, y
1119 >>> factor_terms(Sum(x + y, (x, 1, 3)))
1120 y*Sum(1, (x, 1, 3)) + Sum(x, (x, 1, 3))
1121 >>> factor_terms(Sum(x*y, (x, 1, 3)))
1122 y*Sum(x, (x, 1, 3))
1124 Notes
1125 =====
1127 If a function in the summand or integrand is replaced
1128 with a symbol, then this simplification should not be
1129 done or else an incorrect result will be obtained when
1130 the symbol is replaced with an expression that depends
1131 on the variables of summation/integration:
1133 >>> eq = Sum(y, (x, 1, 3))
1134 >>> factor_terms(eq).subs(y, x).doit()
1135 3*x
1136 >>> eq.subs(y, x).doit()
1137 6
1138 """
1139 result = expr.function
1140 if result == 0:
1141 return S.Zero
1142 limits = expr.limits
1144 # get the wrt variables
1145 wrt = {i.args[0] for i in limits}
1147 # factor out any common terms that are independent of wrt
1148 f = factor_terms(result, **kwargs)
1149 i, d = f.as_independent(*wrt)
1150 if isinstance(f, Add):
1151 return i * expr.func(1, *limits) + expr.func(d, *limits)
1152 else:
1153 return i * expr.func(d, *limits)
1156def factor_terms(expr, radical=False, clear=False, fraction=False, sign=True):
1157 """Remove common factors from terms in all arguments without
1158 changing the underlying structure of the expr. No expansion or
1159 simplification (and no processing of non-commutatives) is performed.
1161 Parameters
1162 ==========
1164 radical: bool, optional
1165 If radical=True then a radical common to all terms will be factored
1166 out of any Add sub-expressions of the expr.
1168 clear : bool, optional
1169 If clear=False (default) then coefficients will not be separated
1170 from a single Add if they can be distributed to leave one or more
1171 terms with integer coefficients.
1173 fraction : bool, optional
1174 If fraction=True (default is False) then a common denominator will be
1175 constructed for the expression.
1177 sign : bool, optional
1178 If sign=True (default) then even if the only factor in common is a -1,
1179 it will be factored out of the expression.
1181 Examples
1182 ========
1184 >>> from sympy import factor_terms, Symbol
1185 >>> from sympy.abc import x, y
1186 >>> factor_terms(x + x*(2 + 4*y)**3)
1187 x*(8*(2*y + 1)**3 + 1)
1188 >>> A = Symbol('A', commutative=False)
1189 >>> factor_terms(x*A + x*A + x*y*A)
1190 x*(y*A + 2*A)
1192 When ``clear`` is False, a rational will only be factored out of an
1193 Add expression if all terms of the Add have coefficients that are
1194 fractions:
1196 >>> factor_terms(x/2 + 1, clear=False)
1197 x/2 + 1
1198 >>> factor_terms(x/2 + 1, clear=True)
1199 (x + 2)/2
1201 If a -1 is all that can be factored out, to *not* factor it out, the
1202 flag ``sign`` must be False:
1204 >>> factor_terms(-x - y)
1205 -(x + y)
1206 >>> factor_terms(-x - y, sign=False)
1207 -x - y
1208 >>> factor_terms(-2*x - 2*y, sign=False)
1209 -2*(x + y)
1211 See Also
1212 ========
1214 gcd_terms, sympy.polys.polytools.terms_gcd
1216 """
1217 def do(expr):
1218 from sympy.concrete.summations import Sum
1219 from sympy.integrals.integrals import Integral
1220 is_iterable = iterable(expr)
1222 if not isinstance(expr, Basic) or expr.is_Atom:
1223 if is_iterable:
1224 return type(expr)([do(i) for i in expr])
1225 return expr
1227 if expr.is_Pow or expr.is_Function or \
1228 is_iterable or not hasattr(expr, 'args_cnc'):
1229 args = expr.args
1230 newargs = tuple([do(i) for i in args])
1231 if newargs == args:
1232 return expr
1233 return expr.func(*newargs)
1235 if isinstance(expr, (Sum, Integral)):
1236 return _factor_sum_int(expr,
1237 radical=radical, clear=clear,
1238 fraction=fraction, sign=sign)
1240 cont, p = expr.as_content_primitive(radical=radical, clear=clear)
1241 if p.is_Add:
1242 list_args = [do(a) for a in Add.make_args(p)]
1243 # get a common negative (if there) which gcd_terms does not remove
1244 if not any(a.as_coeff_Mul()[0].extract_multiplicatively(-1) is None
1245 for a in list_args):
1246 cont = -cont
1247 list_args = [-a for a in list_args]
1248 # watch out for exp(-(x+2)) which gcd_terms will change to exp(-x-2)
1249 special = {}
1250 for i, a in enumerate(list_args):
1251 b, e = a.as_base_exp()
1252 if e.is_Mul and e != Mul(*e.args):
1253 list_args[i] = Dummy()
1254 special[list_args[i]] = a
1255 # rebuild p not worrying about the order which gcd_terms will fix
1256 p = Add._from_args(list_args)
1257 p = gcd_terms(p,
1258 isprimitive=True,
1259 clear=clear,
1260 fraction=fraction).xreplace(special)
1261 elif p.args:
1262 p = p.func(
1263 *[do(a) for a in p.args])
1264 rv = _keep_coeff(cont, p, clear=clear, sign=sign)
1265 return rv
1266 expr = sympify(expr)
1267 return do(expr)
1270def _mask_nc(eq, name=None):
1271 """
1272 Return ``eq`` with non-commutative objects replaced with Dummy
1273 symbols. A dictionary that can be used to restore the original
1274 values is returned: if it is None, the expression is noncommutative
1275 and cannot be made commutative. The third value returned is a list
1276 of any non-commutative symbols that appear in the returned equation.
1278 Explanation
1279 ===========
1281 All non-commutative objects other than Symbols are replaced with
1282 a non-commutative Symbol. Identical objects will be identified
1283 by identical symbols.
1285 If there is only 1 non-commutative object in an expression it will
1286 be replaced with a commutative symbol. Otherwise, the non-commutative
1287 entities are retained and the calling routine should handle
1288 replacements in this case since some care must be taken to keep
1289 track of the ordering of symbols when they occur within Muls.
1291 Parameters
1292 ==========
1294 name : str
1295 ``name``, if given, is the name that will be used with numbered Dummy
1296 variables that will replace the non-commutative objects and is mainly
1297 used for doctesting purposes.
1299 Examples
1300 ========
1302 >>> from sympy.physics.secondquant import Commutator, NO, F, Fd
1303 >>> from sympy import symbols
1304 >>> from sympy.core.exprtools import _mask_nc
1305 >>> from sympy.abc import x, y
1306 >>> A, B, C = symbols('A,B,C', commutative=False)
1308 One nc-symbol:
1310 >>> _mask_nc(A**2 - x**2, 'd')
1311 (_d0**2 - x**2, {_d0: A}, [])
1313 Multiple nc-symbols:
1315 >>> _mask_nc(A**2 - B**2, 'd')
1316 (A**2 - B**2, {}, [A, B])
1318 An nc-object with nc-symbols but no others outside of it:
1320 >>> _mask_nc(1 + x*Commutator(A, B), 'd')
1321 (_d0*x + 1, {_d0: Commutator(A, B)}, [])
1322 >>> _mask_nc(NO(Fd(x)*F(y)), 'd')
1323 (_d0, {_d0: NO(CreateFermion(x)*AnnihilateFermion(y))}, [])
1325 Multiple nc-objects:
1327 >>> eq = x*Commutator(A, B) + x*Commutator(A, C)*Commutator(A, B)
1328 >>> _mask_nc(eq, 'd')
1329 (x*_d0 + x*_d1*_d0, {_d0: Commutator(A, B), _d1: Commutator(A, C)}, [_d0, _d1])
1331 Multiple nc-objects and nc-symbols:
1333 >>> eq = A*Commutator(A, B) + B*Commutator(A, C)
1334 >>> _mask_nc(eq, 'd')
1335 (A*_d0 + B*_d1, {_d0: Commutator(A, B), _d1: Commutator(A, C)}, [_d0, _d1, A, B])
1337 """
1338 name = name or 'mask'
1339 # Make Dummy() append sequential numbers to the name
1341 def numbered_names():
1342 i = 0
1343 while True:
1344 yield name + str(i)
1345 i += 1
1347 names = numbered_names()
1349 def Dummy(*args, **kwargs):
1350 from .symbol import Dummy
1351 return Dummy(next(names), *args, **kwargs)
1353 expr = eq
1354 if expr.is_commutative:
1355 return eq, {}, []
1357 # identify nc-objects; symbols and other
1358 rep = []
1359 nc_obj = set()
1360 nc_syms = set()
1361 pot = preorder_traversal(expr, keys=default_sort_key)
1362 for i, a in enumerate(pot):
1363 if any(a == r[0] for r in rep):
1364 pot.skip()
1365 elif not a.is_commutative:
1366 if a.is_symbol:
1367 nc_syms.add(a)
1368 pot.skip()
1369 elif not (a.is_Add or a.is_Mul or a.is_Pow):
1370 nc_obj.add(a)
1371 pot.skip()
1373 # If there is only one nc symbol or object, it can be factored regularly
1374 # but polys is going to complain, so replace it with a Dummy.
1375 if len(nc_obj) == 1 and not nc_syms:
1376 rep.append((nc_obj.pop(), Dummy()))
1377 elif len(nc_syms) == 1 and not nc_obj:
1378 rep.append((nc_syms.pop(), Dummy()))
1380 # Any remaining nc-objects will be replaced with an nc-Dummy and
1381 # identified as an nc-Symbol to watch out for
1382 nc_obj = sorted(nc_obj, key=default_sort_key)
1383 for n in nc_obj:
1384 nc = Dummy(commutative=False)
1385 rep.append((n, nc))
1386 nc_syms.add(nc)
1387 expr = expr.subs(rep)
1389 nc_syms = list(nc_syms)
1390 nc_syms.sort(key=default_sort_key)
1391 return expr, {v: k for k, v in rep}, nc_syms
1394def factor_nc(expr):
1395 """Return the factored form of ``expr`` while handling non-commutative
1396 expressions.
1398 Examples
1399 ========
1401 >>> from sympy import factor_nc, Symbol
1402 >>> from sympy.abc import x
1403 >>> A = Symbol('A', commutative=False)
1404 >>> B = Symbol('B', commutative=False)
1405 >>> factor_nc((x**2 + 2*A*x + A**2).expand())
1406 (x + A)**2
1407 >>> factor_nc(((x + A)*(x + B)).expand())
1408 (x + A)*(x + B)
1409 """
1410 expr = sympify(expr)
1411 if not isinstance(expr, Expr) or not expr.args:
1412 return expr
1413 if not expr.is_Add:
1414 return expr.func(*[factor_nc(a) for a in expr.args])
1415 expr = expr.func(*[expand_power_exp(i) for i in expr.args])
1417 from sympy.polys.polytools import gcd, factor
1418 expr, rep, nc_symbols = _mask_nc(expr)
1420 if rep:
1421 return factor(expr).subs(rep)
1422 else:
1423 args = [a.args_cnc() for a in Add.make_args(expr)]
1424 c = g = l = r = S.One
1425 hit = False
1426 # find any commutative gcd term
1427 for i, a in enumerate(args):
1428 if i == 0:
1429 c = Mul._from_args(a[0])
1430 elif a[0]:
1431 c = gcd(c, Mul._from_args(a[0]))
1432 else:
1433 c = S.One
1434 if c is not S.One:
1435 hit = True
1436 c, g = c.as_coeff_Mul()
1437 if g is not S.One:
1438 for i, (cc, _) in enumerate(args):
1439 cc = list(Mul.make_args(Mul._from_args(list(cc))/g))
1440 args[i][0] = cc
1441 for i, (cc, _) in enumerate(args):
1442 if cc:
1443 cc[0] = cc[0]/c
1444 else:
1445 cc = [1/c]
1446 args[i][0] = cc
1447 # find any noncommutative common prefix
1448 for i, a in enumerate(args):
1449 if i == 0:
1450 n = a[1][:]
1451 else:
1452 n = common_prefix(n, a[1])
1453 if not n:
1454 # is there a power that can be extracted?
1455 if not args[0][1]:
1456 break
1457 b, e = args[0][1][0].as_base_exp()
1458 ok = False
1459 if e.is_Integer:
1460 for t in args:
1461 if not t[1]:
1462 break
1463 bt, et = t[1][0].as_base_exp()
1464 if et.is_Integer and bt == b:
1465 e = min(e, et)
1466 else:
1467 break
1468 else:
1469 ok = hit = True
1470 l = b**e
1471 il = b**-e
1472 for _ in args:
1473 _[1][0] = il*_[1][0]
1474 break
1475 if not ok:
1476 break
1477 else:
1478 hit = True
1479 lenn = len(n)
1480 l = Mul(*n)
1481 for _ in args:
1482 _[1] = _[1][lenn:]
1483 # find any noncommutative common suffix
1484 for i, a in enumerate(args):
1485 if i == 0:
1486 n = a[1][:]
1487 else:
1488 n = common_suffix(n, a[1])
1489 if not n:
1490 # is there a power that can be extracted?
1491 if not args[0][1]:
1492 break
1493 b, e = args[0][1][-1].as_base_exp()
1494 ok = False
1495 if e.is_Integer:
1496 for t in args:
1497 if not t[1]:
1498 break
1499 bt, et = t[1][-1].as_base_exp()
1500 if et.is_Integer and bt == b:
1501 e = min(e, et)
1502 else:
1503 break
1504 else:
1505 ok = hit = True
1506 r = b**e
1507 il = b**-e
1508 for _ in args:
1509 _[1][-1] = _[1][-1]*il
1510 break
1511 if not ok:
1512 break
1513 else:
1514 hit = True
1515 lenn = len(n)
1516 r = Mul(*n)
1517 for _ in args:
1518 _[1] = _[1][:len(_[1]) - lenn]
1519 if hit:
1520 mid = Add(*[Mul(*cc)*Mul(*nc) for cc, nc in args])
1521 else:
1522 mid = expr
1524 from sympy.simplify.powsimp import powsimp
1526 # sort the symbols so the Dummys would appear in the same
1527 # order as the original symbols, otherwise you may introduce
1528 # a factor of -1, e.g. A**2 - B**2) -- {A:y, B:x} --> y**2 - x**2
1529 # and the former factors into two terms, (A - B)*(A + B) while the
1530 # latter factors into 3 terms, (-1)*(x - y)*(x + y)
1531 rep1 = [(n, Dummy()) for n in sorted(nc_symbols, key=default_sort_key)]
1532 unrep1 = [(v, k) for k, v in rep1]
1533 unrep1.reverse()
1534 new_mid, r2, _ = _mask_nc(mid.subs(rep1))
1535 new_mid = powsimp(factor(new_mid))
1537 new_mid = new_mid.subs(r2).subs(unrep1)
1539 if new_mid.is_Pow:
1540 return _keep_coeff(c, g*l*new_mid*r)
1542 if new_mid.is_Mul:
1543 def _pemexpand(expr):
1544 "Expand with the minimal set of hints necessary to check the result."
1545 return expr.expand(deep=True, mul=True, power_exp=True,
1546 power_base=False, basic=False, multinomial=True, log=False)
1547 # XXX TODO there should be a way to inspect what order the terms
1548 # must be in and just select the plausible ordering without
1549 # checking permutations
1550 cfac = []
1551 ncfac = []
1552 for f in new_mid.args:
1553 if f.is_commutative:
1554 cfac.append(f)
1555 else:
1556 b, e = f.as_base_exp()
1557 if e.is_Integer:
1558 ncfac.extend([b]*e)
1559 else:
1560 ncfac.append(f)
1561 pre_mid = g*Mul(*cfac)*l
1562 target = _pemexpand(expr/c)
1563 for s in variations(ncfac, len(ncfac)):
1564 ok = pre_mid*Mul(*s)*r
1565 if _pemexpand(ok) == target:
1566 return _keep_coeff(c, ok)
1568 # mid was an Add that didn't factor successfully
1569 return _keep_coeff(c, g*l*mid*r)