| 1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573 |
- """Tools for manipulating of large commutative expressions. """
- from __future__ import annotations
- from .add import Add
- from .mul import Mul, _keep_coeff
- from .power import Pow
- from .basic import Basic
- from .expr import Expr
- from .function import expand_power_exp
- from .sympify import sympify
- from .numbers import Rational, Integer, Number, I, equal_valued
- from .singleton import S
- from .sorting import default_sort_key, ordered
- from .symbol import Dummy
- from .traversal import preorder_traversal
- from .coreerrors import NonCommutativeExpression
- from .containers import Tuple, Dict
- from sympy.external.gmpy import SYMPY_INTS
- from sympy.utilities.iterables import (common_prefix, common_suffix,
- variations, iterable, is_sequence)
- from collections import defaultdict
- _eps = Dummy(positive=True)
- def _isnumber(i):
- return isinstance(i, (SYMPY_INTS, float)) or i.is_Number
- def _monotonic_sign(self):
- """Return the value closest to 0 that ``self`` may have if all symbols
- are signed and the result is uniformly the same sign for all values of symbols.
- If a symbol is only signed but not known to be an
- integer or the result is 0 then a symbol representative of the sign of self
- will be returned. Otherwise, None is returned if a) the sign could be positive
- or negative or b) self is not in one of the following forms:
- - L(x, y, ...) + A: a function linear in all symbols x, y, ... with an
- additive constant; if A is zero then the function can be a monomial whose
- sign is monotonic over the range of the variables, e.g. (x + 1)**3 if x is
- nonnegative.
- - A/L(x, y, ...) + B: the inverse of a function linear in all symbols x, y, ...
- that does not have a sign change from positive to negative for any set
- of values for the variables.
- - M(x, y, ...) + A: a monomial M whose factors are all signed and a constant, A.
- - A/M(x, y, ...) + B: the inverse of a monomial and constants A and B.
- - P(x): a univariate polynomial
- Examples
- ========
- >>> from sympy.core.exprtools import _monotonic_sign as F
- >>> from sympy import Dummy
- >>> nn = Dummy(integer=True, nonnegative=True)
- >>> p = Dummy(integer=True, positive=True)
- >>> p2 = Dummy(integer=True, positive=True)
- >>> F(nn + 1)
- 1
- >>> F(p - 1)
- _nneg
- >>> F(nn*p + 1)
- 1
- >>> F(p2*p + 1)
- 2
- >>> F(nn - 1) # could be negative, zero or positive
- """
- if not self.is_extended_real:
- return
- if (-self).is_Symbol:
- rv = _monotonic_sign(-self)
- return rv if rv is None else -rv
- if not self.is_Add and self.as_numer_denom()[1].is_number:
- s = self
- if s.is_prime:
- if s.is_odd:
- return Integer(3)
- else:
- return Integer(2)
- elif s.is_composite:
- if s.is_odd:
- return Integer(9)
- else:
- return Integer(4)
- elif s.is_positive:
- if s.is_even:
- if s.is_prime is False:
- return Integer(4)
- else:
- return Integer(2)
- elif s.is_integer:
- return S.One
- else:
- return _eps
- elif s.is_extended_negative:
- if s.is_even:
- return Integer(-2)
- elif s.is_integer:
- return S.NegativeOne
- else:
- return -_eps
- if s.is_zero or s.is_extended_nonpositive or s.is_extended_nonnegative:
- return S.Zero
- return None
- # univariate polynomial
- free = self.free_symbols
- if len(free) == 1:
- if self.is_polynomial():
- from sympy.polys.polytools import real_roots
- from sympy.polys.polyroots import roots
- from sympy.polys.polyerrors import PolynomialError
- x = free.pop()
- x0 = _monotonic_sign(x)
- if x0 in (_eps, -_eps):
- x0 = S.Zero
- if x0 is not None:
- d = self.diff(x)
- if d.is_number:
- currentroots = []
- else:
- try:
- currentroots = real_roots(d)
- except (PolynomialError, NotImplementedError):
- currentroots = [r for r in roots(d, x) if r.is_extended_real]
- y = self.subs(x, x0)
- if x.is_nonnegative and all(
- (r - x0).is_nonpositive for r in currentroots):
- if y.is_nonnegative and d.is_positive:
- if y:
- return y if y.is_positive else Dummy('pos', positive=True)
- else:
- return Dummy('nneg', nonnegative=True)
- if y.is_nonpositive and d.is_negative:
- if y:
- return y if y.is_negative else Dummy('neg', negative=True)
- else:
- return Dummy('npos', nonpositive=True)
- elif x.is_nonpositive and all(
- (r - x0).is_nonnegative for r in currentroots):
- if y.is_nonnegative and d.is_negative:
- if y:
- return Dummy('pos', positive=True)
- else:
- return Dummy('nneg', nonnegative=True)
- if y.is_nonpositive and d.is_positive:
- if y:
- return Dummy('neg', negative=True)
- else:
- return Dummy('npos', nonpositive=True)
- else:
- n, d = self.as_numer_denom()
- den = None
- if n.is_number:
- den = _monotonic_sign(d)
- elif not d.is_number:
- if _monotonic_sign(n) is not None:
- den = _monotonic_sign(d)
- if den is not None and (den.is_positive or den.is_negative):
- v = n*den
- if v.is_positive:
- return Dummy('pos', positive=True)
- elif v.is_nonnegative:
- return Dummy('nneg', nonnegative=True)
- elif v.is_negative:
- return Dummy('neg', negative=True)
- elif v.is_nonpositive:
- return Dummy('npos', nonpositive=True)
- return None
- # multivariate
- c, a = self.as_coeff_Add()
- v = None
- if not a.is_polynomial():
- # F/A or A/F where A is a number and F is a signed, rational monomial
- n, d = a.as_numer_denom()
- if not (n.is_number or d.is_number):
- return
- if (
- a.is_Mul or a.is_Pow) and \
- a.is_rational and \
- all(p.exp.is_Integer for p in a.atoms(Pow) if p.is_Pow) and \
- (a.is_positive or a.is_negative):
- v = S.One
- for ai in Mul.make_args(a):
- if ai.is_number:
- v *= ai
- continue
- reps = {}
- for x in ai.free_symbols:
- reps[x] = _monotonic_sign(x)
- if reps[x] is None:
- return
- v *= ai.subs(reps)
- elif c:
- # signed linear expression
- if not any(p for p in a.atoms(Pow) if not p.is_number) and (a.is_nonpositive or a.is_nonnegative):
- free = list(a.free_symbols)
- p = {}
- for i in free:
- v = _monotonic_sign(i)
- if v is None:
- return
- p[i] = v or (_eps if i.is_nonnegative else -_eps)
- v = a.xreplace(p)
- if v is not None:
- rv = v + c
- if v.is_nonnegative and rv.is_positive:
- return rv.subs(_eps, 0)
- if v.is_nonpositive and rv.is_negative:
- return rv.subs(_eps, 0)
- def decompose_power(expr: Expr) -> tuple[Expr, int]:
- """
- Decompose power into symbolic base and integer exponent.
- Examples
- ========
- >>> from sympy.core.exprtools import decompose_power
- >>> from sympy.abc import x, y
- >>> from sympy import exp
- >>> decompose_power(x)
- (x, 1)
- >>> decompose_power(x**2)
- (x, 2)
- >>> decompose_power(exp(2*y/3))
- (exp(y/3), 2)
- """
- base, exp = expr.as_base_exp()
- if exp.is_Number:
- if exp.is_Rational:
- if not exp.is_Integer:
- base = Pow(base, Rational(1, exp.q)) # type: ignore
- e = exp.p # type: ignore
- else:
- base, e = expr, 1
- else:
- exp, tail = exp.as_coeff_Mul(rational=True)
- if exp is S.NegativeOne:
- base, e = Pow(base, tail), -1
- elif exp is not S.One:
- # todo: after dropping python 3.7 support, use overload and Literal
- # in as_coeff_Mul to make exp Rational, and remove these 2 ignores
- tail = _keep_coeff(Rational(1, exp.q), tail) # type: ignore
- base, e = Pow(base, tail), exp.p # type: ignore
- else:
- base, e = expr, 1
- return base, e
- def decompose_power_rat(expr: Expr) -> tuple[Expr, Rational]:
- """
- Decompose power into symbolic base and rational exponent;
- if the exponent is not a Rational, then separate only the
- integer coefficient.
- Examples
- ========
- >>> from sympy.core.exprtools import decompose_power_rat
- >>> from sympy.abc import x
- >>> from sympy import sqrt, exp
- >>> decompose_power_rat(sqrt(x))
- (x, 1/2)
- >>> decompose_power_rat(exp(-3*x/2))
- (exp(x/2), -3)
- """
- base, exp = expr.as_base_exp()
- if not exp.is_Rational:
- base, exp_i = decompose_power(expr)
- exp = Integer(exp_i)
- return base, exp # type: ignore
- class Factors:
- """Efficient representation of ``f_1*f_2*...*f_n``."""
- __slots__ = ('factors', 'gens')
- def __init__(self, factors=None): # Factors
- """Initialize Factors from dict or expr.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x
- >>> from sympy import I
- >>> e = 2*x**3
- >>> Factors(e)
- Factors({2: 1, x: 3})
- >>> Factors(e.as_powers_dict())
- Factors({2: 1, x: 3})
- >>> f = _
- >>> f.factors # underlying dictionary
- {2: 1, x: 3}
- >>> f.gens # base of each factor
- frozenset({2, x})
- >>> Factors(0)
- Factors({0: 1})
- >>> Factors(I)
- Factors({I: 1})
- Notes
- =====
- Although a dictionary can be passed, only minimal checking is
- performed: powers of -1 and I are made canonical.
- """
- if isinstance(factors, (SYMPY_INTS, float)):
- factors = S(factors)
- if isinstance(factors, Factors):
- factors = factors.factors.copy()
- elif factors in (None, S.One):
- factors = {}
- elif factors is S.Zero or factors == 0:
- factors = {S.Zero: S.One}
- elif isinstance(factors, Number):
- n = factors
- factors = {}
- if n < 0:
- factors[S.NegativeOne] = S.One
- n = -n
- if n is not S.One:
- if n.is_Float or n.is_Integer or n is S.Infinity:
- factors[n] = S.One
- elif n.is_Rational:
- # since we're processing Numbers, the denominator is
- # stored with a negative exponent; all other factors
- # are left .
- if n.p != 1:
- factors[Integer(n.p)] = S.One
- factors[Integer(n.q)] = S.NegativeOne
- else:
- raise ValueError('Expected Float|Rational|Integer, not %s' % n)
- elif isinstance(factors, Basic) and not factors.args:
- factors = {factors: S.One}
- elif isinstance(factors, Expr):
- c, nc = factors.args_cnc()
- i = c.count(I)
- for _ in range(i):
- c.remove(I)
- factors = dict(Mul._from_args(c).as_powers_dict())
- # Handle all rational Coefficients
- for f in list(factors.keys()):
- if isinstance(f, Rational) and not isinstance(f, Integer):
- p, q = Integer(f.p), Integer(f.q)
- factors[p] = (factors[p] if p in factors else S.Zero) + factors[f]
- factors[q] = (factors[q] if q in factors else S.Zero) - factors[f]
- factors.pop(f)
- if i:
- factors[I] = factors.get(I, S.Zero) + i
- if nc:
- factors[Mul(*nc, evaluate=False)] = S.One
- else:
- factors = factors.copy() # /!\ should be dict-like
- # tidy up -/+1 and I exponents if Rational
- handle = [k for k in factors if k is I or k in (-1, 1)]
- if handle:
- i1 = S.One
- for k in handle:
- if not _isnumber(factors[k]):
- continue
- i1 *= k**factors.pop(k)
- if i1 is not S.One:
- for a in i1.args if i1.is_Mul else [i1]: # at worst, -1.0*I*(-1)**e
- if a is S.NegativeOne:
- factors[a] = S.One
- elif a is I:
- factors[I] = S.One
- elif a.is_Pow:
- factors[a.base] = factors.get(a.base, S.Zero) + a.exp
- elif equal_valued(a, 1):
- factors[a] = S.One
- elif equal_valued(a, -1):
- factors[-a] = S.One
- factors[S.NegativeOne] = S.One
- else:
- raise ValueError('unexpected factor in i1: %s' % a)
- self.factors = factors
- keys = getattr(factors, 'keys', None)
- if keys is None:
- raise TypeError('expecting Expr or dictionary')
- self.gens = frozenset(keys())
- def __hash__(self): # Factors
- keys = tuple(ordered(self.factors.keys()))
- values = [self.factors[k] for k in keys]
- return hash((keys, values))
- def __repr__(self): # Factors
- return "Factors({%s})" % ', '.join(
- ['%s: %s' % (k, v) for k, v in ordered(self.factors.items())])
- @property
- def is_zero(self): # Factors
- """
- >>> from sympy.core.exprtools import Factors
- >>> Factors(0).is_zero
- True
- """
- f = self.factors
- return len(f) == 1 and S.Zero in f
- @property
- def is_one(self): # Factors
- """
- >>> from sympy.core.exprtools import Factors
- >>> Factors(1).is_one
- True
- """
- return not self.factors
- def as_expr(self): # Factors
- """Return the underlying expression.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y
- >>> Factors((x*y**2).as_powers_dict()).as_expr()
- x*y**2
- """
- args = []
- for factor, exp in self.factors.items():
- if exp != 1:
- if isinstance(exp, Integer):
- b, e = factor.as_base_exp()
- e = _keep_coeff(exp, e)
- args.append(b**e)
- else:
- args.append(factor**exp)
- else:
- args.append(factor)
- return Mul(*args)
- def mul(self, other): # Factors
- """Return Factors of ``self * other``.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y, z
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> b = Factors((x*y/z).as_powers_dict())
- >>> a.mul(b)
- Factors({x: 2, y: 3, z: -1})
- >>> a*b
- Factors({x: 2, y: 3, z: -1})
- """
- if not isinstance(other, Factors):
- other = Factors(other)
- if any(f.is_zero for f in (self, other)):
- return Factors(S.Zero)
- factors = dict(self.factors)
- for factor, exp in other.factors.items():
- if factor in factors:
- exp = factors[factor] + exp
- if not exp:
- del factors[factor]
- continue
- factors[factor] = exp
- return Factors(factors)
- def normal(self, other):
- """Return ``self`` and ``other`` with ``gcd`` removed from each.
- The only differences between this and method ``div`` is that this
- is 1) optimized for the case when there are few factors in common and
- 2) this does not raise an error if ``other`` is zero.
- See Also
- ========
- div
- """
- if not isinstance(other, Factors):
- other = Factors(other)
- if other.is_zero:
- return (Factors(), Factors(S.Zero))
- if self.is_zero:
- return (Factors(S.Zero), Factors())
- self_factors = dict(self.factors)
- other_factors = dict(other.factors)
- for factor, self_exp in self.factors.items():
- try:
- other_exp = other.factors[factor]
- except KeyError:
- continue
- exp = self_exp - other_exp
- if not exp:
- del self_factors[factor]
- del other_factors[factor]
- elif _isnumber(exp):
- if exp > 0:
- self_factors[factor] = exp
- del other_factors[factor]
- else:
- del self_factors[factor]
- other_factors[factor] = -exp
- else:
- r = self_exp.extract_additively(other_exp)
- if r is not None:
- if r:
- self_factors[factor] = r
- del other_factors[factor]
- else: # should be handled already
- del self_factors[factor]
- del other_factors[factor]
- else:
- sc, sa = self_exp.as_coeff_Add()
- if sc:
- oc, oa = other_exp.as_coeff_Add()
- diff = sc - oc
- if diff > 0:
- self_factors[factor] -= oc
- other_exp = oa
- elif diff < 0:
- self_factors[factor] -= sc
- other_factors[factor] -= sc
- other_exp = oa - diff
- else:
- self_factors[factor] = sa
- other_exp = oa
- if other_exp:
- other_factors[factor] = other_exp
- else:
- del other_factors[factor]
- return Factors(self_factors), Factors(other_factors)
- def div(self, other): # Factors
- """Return ``self`` and ``other`` with ``gcd`` removed from each.
- This is optimized for the case when there are many factors in common.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y, z
- >>> from sympy import S
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> a.div(a)
- (Factors({}), Factors({}))
- >>> a.div(x*z)
- (Factors({y: 2}), Factors({z: 1}))
- The ``/`` operator only gives ``quo``:
- >>> a/x
- Factors({y: 2})
- Factors treats its factors as though they are all in the numerator, so
- if you violate this assumption the results will be correct but will
- not strictly correspond to the numerator and denominator of the ratio:
- >>> a.div(x/z)
- (Factors({y: 2}), Factors({z: -1}))
- Factors is also naive about bases: it does not attempt any denesting
- of Rational-base terms, for example the following does not become
- 2**(2*x)/2.
- >>> Factors(2**(2*x + 2)).div(S(8))
- (Factors({2: 2*x + 2}), Factors({8: 1}))
- factor_terms can clean up such Rational-bases powers:
- >>> from sympy import factor_terms
- >>> n, d = Factors(2**(2*x + 2)).div(S(8))
- >>> n.as_expr()/d.as_expr()
- 2**(2*x + 2)/8
- >>> factor_terms(_)
- 2**(2*x)/2
- """
- quo, rem = dict(self.factors), {}
- if not isinstance(other, Factors):
- other = Factors(other)
- if other.is_zero:
- raise ZeroDivisionError
- if self.is_zero:
- return (Factors(S.Zero), Factors())
- for factor, exp in other.factors.items():
- if factor in quo:
- d = quo[factor] - exp
- if _isnumber(d):
- if d <= 0:
- del quo[factor]
- if d >= 0:
- if d:
- quo[factor] = d
- continue
- exp = -d
- else:
- r = quo[factor].extract_additively(exp)
- if r is not None:
- if r:
- quo[factor] = r
- else: # should be handled already
- del quo[factor]
- else:
- other_exp = exp
- sc, sa = quo[factor].as_coeff_Add()
- if sc:
- oc, oa = other_exp.as_coeff_Add()
- diff = sc - oc
- if diff > 0:
- quo[factor] -= oc
- other_exp = oa
- elif diff < 0:
- quo[factor] -= sc
- other_exp = oa - diff
- else:
- quo[factor] = sa
- other_exp = oa
- if other_exp:
- rem[factor] = other_exp
- else:
- assert factor not in rem
- continue
- rem[factor] = exp
- return Factors(quo), Factors(rem)
- def quo(self, other): # Factors
- """Return numerator Factor of ``self / other``.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y, z
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> b = Factors((x*y/z).as_powers_dict())
- >>> a.quo(b) # same as a/b
- Factors({y: 1})
- """
- return self.div(other)[0]
- def rem(self, other): # Factors
- """Return denominator Factors of ``self / other``.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y, z
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> b = Factors((x*y/z).as_powers_dict())
- >>> a.rem(b)
- Factors({z: -1})
- >>> a.rem(a)
- Factors({})
- """
- return self.div(other)[1]
- def pow(self, other): # Factors
- """Return self raised to a non-negative integer power.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> a**2
- Factors({x: 2, y: 4})
- """
- if isinstance(other, Factors):
- other = other.as_expr()
- if other.is_Integer:
- other = int(other)
- if isinstance(other, SYMPY_INTS) and other >= 0:
- factors = {}
- if other:
- for factor, exp in self.factors.items():
- factors[factor] = exp*other
- return Factors(factors)
- else:
- raise ValueError("expected non-negative integer, got %s" % other)
- def gcd(self, other): # Factors
- """Return Factors of ``gcd(self, other)``. The keys are
- the intersection of factors with the minimum exponent for
- each factor.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y, z
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> b = Factors((x*y/z).as_powers_dict())
- >>> a.gcd(b)
- Factors({x: 1, y: 1})
- """
- if not isinstance(other, Factors):
- other = Factors(other)
- if other.is_zero:
- return Factors(self.factors)
- factors = {}
- for factor, exp in self.factors.items():
- factor, exp = sympify(factor), sympify(exp)
- if factor in other.factors:
- lt = (exp - other.factors[factor]).is_negative
- if lt == True:
- factors[factor] = exp
- elif lt == False:
- factors[factor] = other.factors[factor]
- return Factors(factors)
- def lcm(self, other): # Factors
- """Return Factors of ``lcm(self, other)`` which are
- the union of factors with the maximum exponent for
- each factor.
- Examples
- ========
- >>> from sympy.core.exprtools import Factors
- >>> from sympy.abc import x, y, z
- >>> a = Factors((x*y**2).as_powers_dict())
- >>> b = Factors((x*y/z).as_powers_dict())
- >>> a.lcm(b)
- Factors({x: 1, y: 2, z: -1})
- """
- if not isinstance(other, Factors):
- other = Factors(other)
- if any(f.is_zero for f in (self, other)):
- return Factors(S.Zero)
- factors = dict(self.factors)
- for factor, exp in other.factors.items():
- if factor in factors:
- exp = max(exp, factors[factor])
- factors[factor] = exp
- return Factors(factors)
- def __mul__(self, other): # Factors
- return self.mul(other)
- def __divmod__(self, other): # Factors
- return self.div(other)
- def __truediv__(self, other): # Factors
- return self.quo(other)
- def __mod__(self, other): # Factors
- return self.rem(other)
- def __pow__(self, other): # Factors
- return self.pow(other)
- def __eq__(self, other): # Factors
- if not isinstance(other, Factors):
- other = Factors(other)
- return self.factors == other.factors
- def __ne__(self, other): # Factors
- return not self == other
- class Term:
- """Efficient representation of ``coeff*(numer/denom)``. """
- __slots__ = ('coeff', 'numer', 'denom')
- def __init__(self, term, numer=None, denom=None): # Term
- if numer is None and denom is None:
- if not term.is_commutative:
- raise NonCommutativeExpression(
- 'commutative expression expected')
- coeff, factors = term.as_coeff_mul()
- numer, denom = defaultdict(int), defaultdict(int)
- for factor in factors:
- base, exp = decompose_power(factor)
- if base.is_Add:
- cont, base = base.primitive()
- coeff *= cont**exp
- if exp > 0:
- numer[base] += exp
- else:
- denom[base] += -exp
- numer = Factors(numer)
- denom = Factors(denom)
- else:
- coeff = term
- if numer is None:
- numer = Factors()
- if denom is None:
- denom = Factors()
- self.coeff = coeff
- self.numer = numer
- self.denom = denom
- def __hash__(self): # Term
- return hash((self.coeff, self.numer, self.denom))
- def __repr__(self): # Term
- return "Term(%s, %s, %s)" % (self.coeff, self.numer, self.denom)
- def as_expr(self): # Term
- return self.coeff*(self.numer.as_expr()/self.denom.as_expr())
- def mul(self, other): # Term
- coeff = self.coeff*other.coeff
- numer = self.numer.mul(other.numer)
- denom = self.denom.mul(other.denom)
- numer, denom = numer.normal(denom)
- return Term(coeff, numer, denom)
- def inv(self): # Term
- return Term(1/self.coeff, self.denom, self.numer)
- def quo(self, other): # Term
- return self.mul(other.inv())
- def pow(self, other): # Term
- if other < 0:
- return self.inv().pow(-other)
- else:
- return Term(self.coeff ** other,
- self.numer.pow(other),
- self.denom.pow(other))
- def gcd(self, other): # Term
- return Term(self.coeff.gcd(other.coeff),
- self.numer.gcd(other.numer),
- self.denom.gcd(other.denom))
- def lcm(self, other): # Term
- return Term(self.coeff.lcm(other.coeff),
- self.numer.lcm(other.numer),
- self.denom.lcm(other.denom))
- def __mul__(self, other): # Term
- if isinstance(other, Term):
- return self.mul(other)
- else:
- return NotImplemented
- def __truediv__(self, other): # Term
- if isinstance(other, Term):
- return self.quo(other)
- else:
- return NotImplemented
- def __pow__(self, other): # Term
- if isinstance(other, SYMPY_INTS):
- return self.pow(other)
- else:
- return NotImplemented
- def __eq__(self, other): # Term
- return (self.coeff == other.coeff and
- self.numer == other.numer and
- self.denom == other.denom)
- def __ne__(self, other): # Term
- return not self == other
- def _gcd_terms(terms, isprimitive=False, fraction=True):
- """Helper function for :func:`gcd_terms`.
- Parameters
- ==========
- isprimitive : boolean, optional
- If ``isprimitive`` is True then the call to primitive
- for an Add will be skipped. This is useful when the
- content has already been extracted.
- fraction : boolean, optional
- If ``fraction`` is True then the expression will appear over a common
- denominator, the lcm of all term denominators.
- """
- if isinstance(terms, Basic) and not isinstance(terms, Tuple):
- terms = Add.make_args(terms)
- terms = list(map(Term, [t for t in terms if t]))
- # there is some simplification that may happen if we leave this
- # here rather than duplicate it before the mapping of Term onto
- # the terms
- if len(terms) == 0:
- return S.Zero, S.Zero, S.One
- if len(terms) == 1:
- cont = terms[0].coeff
- numer = terms[0].numer.as_expr()
- denom = terms[0].denom.as_expr()
- else:
- cont = terms[0]
- for term in terms[1:]:
- cont = cont.gcd(term)
- for i, term in enumerate(terms):
- terms[i] = term.quo(cont)
- if fraction:
- denom = terms[0].denom
- for term in terms[1:]:
- denom = denom.lcm(term.denom)
- numers = []
- for term in terms:
- numer = term.numer.mul(denom.quo(term.denom))
- numers.append(term.coeff*numer.as_expr())
- else:
- numers = [t.as_expr() for t in terms]
- denom = Term(S.One).numer
- cont = cont.as_expr()
- numer = Add(*numers)
- denom = denom.as_expr()
- if not isprimitive and numer.is_Add:
- _cont, numer = numer.primitive()
- cont *= _cont
- return cont, numer, denom
- def gcd_terms(terms, isprimitive=False, clear=True, fraction=True):
- """Compute the GCD of ``terms`` and put them together.
- Parameters
- ==========
- terms : Expr
- Can be an expression or a non-Basic sequence of expressions
- which will be handled as though they are terms from a sum.
- isprimitive : bool, optional
- If ``isprimitive`` is True the _gcd_terms will not run the primitive
- method on the terms.
- clear : bool, optional
- It controls the removal of integers from the denominator of an Add
- expression. When True (default), all numerical denominator will be cleared;
- when False the denominators will be cleared only if all terms had numerical
- denominators other than 1.
- fraction : bool, optional
- When True (default), will put the expression over a common
- denominator.
- Examples
- ========
- >>> from sympy import gcd_terms
- >>> from sympy.abc import x, y
- >>> gcd_terms((x + 1)**2*y + (x + 1)*y**2)
- y*(x + 1)*(x + y + 1)
- >>> gcd_terms(x/2 + 1)
- (x + 2)/2
- >>> gcd_terms(x/2 + 1, clear=False)
- x/2 + 1
- >>> gcd_terms(x/2 + y/2, clear=False)
- (x + y)/2
- >>> gcd_terms(x/2 + 1/x)
- (x**2 + 2)/(2*x)
- >>> gcd_terms(x/2 + 1/x, fraction=False)
- (x + 2/x)/2
- >>> gcd_terms(x/2 + 1/x, fraction=False, clear=False)
- x/2 + 1/x
- >>> gcd_terms(x/2/y + 1/x/y)
- (x**2 + 2)/(2*x*y)
- >>> gcd_terms(x/2/y + 1/x/y, clear=False)
- (x**2/2 + 1)/(x*y)
- >>> gcd_terms(x/2/y + 1/x/y, clear=False, fraction=False)
- (x/2 + 1/x)/y
- The ``clear`` flag was ignored in this case because the returned
- expression was a rational expression, not a simple sum.
- See Also
- ========
- factor_terms, sympy.polys.polytools.terms_gcd
- """
- def mask(terms):
- """replace nc portions of each term with a unique Dummy symbols
- and return the replacements to restore them"""
- args = [(a, []) if a.is_commutative else a.args_cnc() for a in terms]
- reps = []
- for i, (c, nc) in enumerate(args):
- if nc:
- nc = Mul(*nc)
- d = Dummy()
- reps.append((d, nc))
- c.append(d)
- args[i] = Mul(*c)
- else:
- args[i] = c
- return args, dict(reps)
- isadd = isinstance(terms, Add)
- addlike = isadd or not isinstance(terms, Basic) and \
- is_sequence(terms, include=set) and \
- not isinstance(terms, Dict)
- if addlike:
- if isadd: # i.e. an Add
- terms = list(terms.args)
- else:
- terms = sympify(terms)
- terms, reps = mask(terms)
- cont, numer, denom = _gcd_terms(terms, isprimitive, fraction)
- numer = numer.xreplace(reps)
- coeff, factors = cont.as_coeff_Mul()
- if not clear:
- c, _coeff = coeff.as_coeff_Mul()
- if not c.is_Integer and not clear and numer.is_Add:
- n, d = c.as_numer_denom()
- _numer = numer/d
- if any(a.as_coeff_Mul()[0].is_Integer
- for a in _numer.args):
- numer = _numer
- coeff = n*_coeff
- return _keep_coeff(coeff, factors*numer/denom, clear=clear)
- if not isinstance(terms, Basic):
- return terms
- if terms.is_Atom:
- return terms
- if terms.is_Mul:
- c, args = terms.as_coeff_mul()
- return _keep_coeff(c, Mul(*[gcd_terms(i, isprimitive, clear, fraction)
- for i in args]), clear=clear)
- def handle(a):
- # don't treat internal args like terms of an Add
- if not isinstance(a, Expr):
- if isinstance(a, Basic):
- if not a.args:
- return a
- return a.func(*[handle(i) for i in a.args])
- return type(a)([handle(i) for i in a])
- return gcd_terms(a, isprimitive, clear, fraction)
- if isinstance(terms, Dict):
- return Dict(*[(k, handle(v)) for k, v in terms.args])
- return terms.func(*[handle(i) for i in terms.args])
- def _factor_sum_int(expr, **kwargs):
- """Return Sum or Integral object with factors that are not
- in the wrt variables removed. In cases where there are additive
- terms in the function of the object that are independent, the
- object will be separated into two objects.
- Examples
- ========
- >>> from sympy import Sum, factor_terms
- >>> from sympy.abc import x, y
- >>> factor_terms(Sum(x + y, (x, 1, 3)))
- y*Sum(1, (x, 1, 3)) + Sum(x, (x, 1, 3))
- >>> factor_terms(Sum(x*y, (x, 1, 3)))
- y*Sum(x, (x, 1, 3))
- Notes
- =====
- If a function in the summand or integrand is replaced
- with a symbol, then this simplification should not be
- done or else an incorrect result will be obtained when
- the symbol is replaced with an expression that depends
- on the variables of summation/integration:
- >>> eq = Sum(y, (x, 1, 3))
- >>> factor_terms(eq).subs(y, x).doit()
- 3*x
- >>> eq.subs(y, x).doit()
- 6
- """
- result = expr.function
- if result == 0:
- return S.Zero
- limits = expr.limits
- # get the wrt variables
- wrt = {i.args[0] for i in limits}
- # factor out any common terms that are independent of wrt
- f = factor_terms(result, **kwargs)
- i, d = f.as_independent(*wrt)
- if isinstance(f, Add):
- return i * expr.func(1, *limits) + expr.func(d, *limits)
- else:
- return i * expr.func(d, *limits)
- def factor_terms(expr: Expr | complex, radical=False, clear=False, fraction=False, sign=True) -> Expr:
- """Remove common factors from terms in all arguments without
- changing the underlying structure of the expr. No expansion or
- simplification (and no processing of non-commutatives) is performed.
- Parameters
- ==========
- radical: bool, optional
- If radical=True then a radical common to all terms will be factored
- out of any Add sub-expressions of the expr.
- clear : bool, optional
- If clear=False (default) then coefficients will not be separated
- from a single Add if they can be distributed to leave one or more
- terms with integer coefficients.
- fraction : bool, optional
- If fraction=True (default is False) then a common denominator will be
- constructed for the expression.
- sign : bool, optional
- If sign=True (default) then even if the only factor in common is a -1,
- it will be factored out of the expression.
- Examples
- ========
- >>> from sympy import factor_terms, Symbol
- >>> from sympy.abc import x, y
- >>> factor_terms(x + x*(2 + 4*y)**3)
- x*(8*(2*y + 1)**3 + 1)
- >>> A = Symbol('A', commutative=False)
- >>> factor_terms(x*A + x*A + x*y*A)
- x*(y*A + 2*A)
- When ``clear`` is False, a rational will only be factored out of an
- Add expression if all terms of the Add have coefficients that are
- fractions:
- >>> factor_terms(x/2 + 1, clear=False)
- x/2 + 1
- >>> factor_terms(x/2 + 1, clear=True)
- (x + 2)/2
- If a -1 is all that can be factored out, to *not* factor it out, the
- flag ``sign`` must be False:
- >>> factor_terms(-x - y)
- -(x + y)
- >>> factor_terms(-x - y, sign=False)
- -x - y
- >>> factor_terms(-2*x - 2*y, sign=False)
- -2*(x + y)
- See Also
- ========
- gcd_terms, sympy.polys.polytools.terms_gcd
- """
- def do(expr):
- from sympy.concrete.summations import Sum
- from sympy.integrals.integrals import Integral
- is_iterable = iterable(expr)
- if not isinstance(expr, Basic) or expr.is_Atom:
- if is_iterable:
- return type(expr)([do(i) for i in expr])
- return expr
- if expr.is_Pow or expr.is_Function or \
- is_iterable or not hasattr(expr, 'args_cnc'):
- args = expr.args
- newargs = tuple([do(i) for i in args])
- if newargs == args:
- return expr
- return expr.func(*newargs)
- if isinstance(expr, (Sum, Integral)):
- return _factor_sum_int(expr,
- radical=radical, clear=clear,
- fraction=fraction, sign=sign)
- cont, p = expr.as_content_primitive(radical=radical, clear=clear)
- if p.is_Add:
- list_args = [do(a) for a in Add.make_args(p)]
- # get a common negative (if there) which gcd_terms does not remove
- if not any(a.as_coeff_Mul()[0].extract_multiplicatively(-1) is None
- for a in list_args):
- cont = -cont
- list_args = [-a for a in list_args]
- # watch out for exp(-(x+2)) which gcd_terms will change to exp(-x-2)
- special = {}
- for i, a in enumerate(list_args):
- b, e = a.as_base_exp()
- if e.is_Mul and e != Mul(*e.args):
- list_args[i] = Dummy()
- special[list_args[i]] = a
- # rebuild p not worrying about the order which gcd_terms will fix
- p = Add._from_args(list_args)
- p = gcd_terms(p,
- isprimitive=True,
- clear=clear,
- fraction=fraction).xreplace(special)
- elif p.args:
- p = p.func(
- *[do(a) for a in p.args])
- rv = _keep_coeff(cont, p, clear=clear, sign=sign)
- return rv
- expr2 = sympify(expr)
- return do(expr2)
- def _mask_nc(eq, name=None):
- """
- Return ``eq`` with non-commutative objects replaced with Dummy
- symbols. A dictionary that can be used to restore the original
- values is returned: if it is None, the expression is noncommutative
- and cannot be made commutative. The third value returned is a list
- of any non-commutative symbols that appear in the returned equation.
- Explanation
- ===========
- All non-commutative objects other than Symbols are replaced with
- a non-commutative Symbol. Identical objects will be identified
- by identical symbols.
- If there is only 1 non-commutative object in an expression it will
- be replaced with a commutative symbol. Otherwise, the non-commutative
- entities are retained and the calling routine should handle
- replacements in this case since some care must be taken to keep
- track of the ordering of symbols when they occur within Muls.
- Parameters
- ==========
- name : str
- ``name``, if given, is the name that will be used with numbered Dummy
- variables that will replace the non-commutative objects and is mainly
- used for doctesting purposes.
- Examples
- ========
- >>> from sympy.physics.secondquant import Commutator, NO, F, Fd
- >>> from sympy import symbols
- >>> from sympy.core.exprtools import _mask_nc
- >>> from sympy.abc import x, y
- >>> A, B, C = symbols('A,B,C', commutative=False)
- One nc-symbol:
- >>> _mask_nc(A**2 - x**2, 'd')
- (_d0**2 - x**2, {_d0: A}, [])
- Multiple nc-symbols:
- >>> _mask_nc(A**2 - B**2, 'd')
- (A**2 - B**2, {}, [A, B])
- An nc-object with nc-symbols but no others outside of it:
- >>> _mask_nc(1 + x*Commutator(A, B), 'd')
- (_d0*x + 1, {_d0: Commutator(A, B)}, [])
- >>> _mask_nc(NO(Fd(x)*F(y)), 'd')
- (_d0, {_d0: NO(CreateFermion(x)*AnnihilateFermion(y))}, [])
- Multiple nc-objects:
- >>> eq = x*Commutator(A, B) + x*Commutator(A, C)*Commutator(A, B)
- >>> _mask_nc(eq, 'd')
- (x*_d0 + x*_d1*_d0, {_d0: Commutator(A, B), _d1: Commutator(A, C)}, [_d0, _d1])
- Multiple nc-objects and nc-symbols:
- >>> eq = A*Commutator(A, B) + B*Commutator(A, C)
- >>> _mask_nc(eq, 'd')
- (A*_d0 + B*_d1, {_d0: Commutator(A, B), _d1: Commutator(A, C)}, [_d0, _d1, A, B])
- """
- name = name or 'mask'
- # Make Dummy() append sequential numbers to the name
- def numbered_names():
- i = 0
- while True:
- yield name + str(i)
- i += 1
- names = numbered_names()
- def Dummy(*args, **kwargs):
- from .symbol import Dummy
- return Dummy(next(names), *args, **kwargs)
- expr = eq
- if expr.is_commutative:
- return eq, {}, []
- # identify nc-objects; symbols and other
- rep = []
- nc_obj = set()
- nc_syms = set()
- pot = preorder_traversal(expr, keys=default_sort_key)
- for a in pot:
- if any(a == r[0] for r in rep):
- pot.skip()
- elif not a.is_commutative:
- if a.is_symbol:
- nc_syms.add(a)
- pot.skip()
- elif not (a.is_Add or a.is_Mul or a.is_Pow):
- nc_obj.add(a)
- pot.skip()
- # If there is only one nc symbol or object, it can be factored regularly
- # but polys is going to complain, so replace it with a Dummy.
- if len(nc_obj) == 1 and not nc_syms:
- rep.append((nc_obj.pop(), Dummy()))
- elif len(nc_syms) == 1 and not nc_obj:
- rep.append((nc_syms.pop(), Dummy()))
- # Any remaining nc-objects will be replaced with an nc-Dummy and
- # identified as an nc-Symbol to watch out for
- nc_obj = sorted(nc_obj, key=default_sort_key)
- for n in nc_obj:
- nc = Dummy(commutative=False)
- rep.append((n, nc))
- nc_syms.add(nc)
- expr = expr.subs(rep)
- nc_syms = list(nc_syms)
- nc_syms.sort(key=default_sort_key)
- return expr, {v: k for k, v in rep}, nc_syms
- def factor_nc(expr):
- """Return the factored form of ``expr`` while handling non-commutative
- expressions.
- Examples
- ========
- >>> from sympy import factor_nc, Symbol
- >>> from sympy.abc import x
- >>> A = Symbol('A', commutative=False)
- >>> B = Symbol('B', commutative=False)
- >>> factor_nc((x**2 + 2*A*x + A**2).expand())
- (x + A)**2
- >>> factor_nc(((x + A)*(x + B)).expand())
- (x + A)*(x + B)
- """
- expr = sympify(expr)
- if not isinstance(expr, Expr) or not expr.args:
- return expr
- if not expr.is_Add:
- return expr.func(*[factor_nc(a) for a in expr.args])
- expr = expr.func(*[expand_power_exp(i) for i in expr.args])
- from sympy.polys.polytools import gcd, factor
- expr, rep, nc_symbols = _mask_nc(expr)
- if rep:
- return factor(expr).subs(rep)
- else:
- args = [a.args_cnc() for a in Add.make_args(expr)]
- c = g = l = r = S.One
- hit = False
- # find any commutative gcd term
- for i, a in enumerate(args):
- if i == 0:
- c = Mul._from_args(a[0])
- elif a[0]:
- c = gcd(c, Mul._from_args(a[0]))
- else:
- c = S.One
- if c is not S.One:
- hit = True
- c, g = c.as_coeff_Mul()
- if g is not S.One:
- for i, (cc, _) in enumerate(args):
- cc = list(Mul.make_args(Mul._from_args(list(cc))/g))
- args[i][0] = cc
- for i, (cc, _) in enumerate(args):
- if cc:
- cc[0] = cc[0]/c
- else:
- cc = [1/c]
- args[i][0] = cc
- # find any noncommutative common prefix
- for i, a in enumerate(args):
- if i == 0:
- n = a[1][:]
- else:
- n = common_prefix(n, a[1])
- if not n:
- # is there a power that can be extracted?
- if not args[0][1]:
- break
- b, e = args[0][1][0].as_base_exp()
- ok = False
- if e.is_Integer:
- for t in args:
- if not t[1]:
- break
- bt, et = t[1][0].as_base_exp()
- if et.is_Integer and bt == b:
- e = min(e, et)
- else:
- break
- else:
- ok = hit = True
- l = b**e
- il = b**-e
- for _ in args:
- _[1][0] = il*_[1][0]
- break
- if not ok:
- break
- else:
- hit = True
- lenn = len(n)
- l = Mul(*n)
- for _ in args:
- _[1] = _[1][lenn:]
- # find any noncommutative common suffix
- for i, a in enumerate(args):
- if i == 0:
- n = a[1][:]
- else:
- n = common_suffix(n, a[1])
- if not n:
- # is there a power that can be extracted?
- if not args[0][1]:
- break
- b, e = args[0][1][-1].as_base_exp()
- ok = False
- if e.is_Integer:
- for t in args:
- if not t[1]:
- break
- bt, et = t[1][-1].as_base_exp()
- if et.is_Integer and bt == b:
- e = min(e, et)
- else:
- break
- else:
- ok = hit = True
- r = b**e
- il = b**-e
- for _ in args:
- _[1][-1] = _[1][-1]*il
- break
- if not ok:
- break
- else:
- hit = True
- lenn = len(n)
- r = Mul(*n)
- for _ in args:
- _[1] = _[1][:len(_[1]) - lenn]
- if hit:
- mid = Add(*[Mul(*cc)*Mul(*nc) for cc, nc in args])
- else:
- mid = expr
- from sympy.simplify.powsimp import powsimp
- # sort the symbols so the Dummys would appear in the same
- # order as the original symbols, otherwise you may introduce
- # a factor of -1, e.g. A**2 - B**2) -- {A:y, B:x} --> y**2 - x**2
- # and the former factors into two terms, (A - B)*(A + B) while the
- # latter factors into 3 terms, (-1)*(x - y)*(x + y)
- rep1 = [(n, Dummy()) for n in sorted(nc_symbols, key=default_sort_key)]
- unrep1 = [(v, k) for k, v in rep1]
- unrep1.reverse()
- new_mid, r2, _ = _mask_nc(mid.subs(rep1))
- new_mid = powsimp(factor(new_mid))
- new_mid = new_mid.subs(r2).subs(unrep1)
- if new_mid.is_Pow:
- return _keep_coeff(c, g*l*new_mid*r)
- if new_mid.is_Mul:
- def _pemexpand(expr):
- "Expand with the minimal set of hints necessary to check the result."
- return expr.expand(deep=True, mul=True, power_exp=True,
- power_base=False, basic=False, multinomial=True, log=False)
- # XXX TODO there should be a way to inspect what order the terms
- # must be in and just select the plausible ordering without
- # checking permutations
- cfac = []
- ncfac = []
- for f in new_mid.args:
- if f.is_commutative:
- cfac.append(f)
- else:
- b, e = f.as_base_exp()
- if e.is_Integer:
- ncfac.extend([b]*e)
- else:
- ncfac.append(f)
- pre_mid = g*Mul(*cfac)*l
- target = _pemexpand(expr/c)
- for s in variations(ncfac, len(ncfac)):
- ok = pre_mid*Mul(*s)*r
- if _pemexpand(ok) == target:
- return _keep_coeff(c, ok)
- # mid was an Add that didn't factor successfully
- return _keep_coeff(c, g*l*mid*r)
|