exprtools.py 50 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573
  1. """Tools for manipulating of large commutative expressions. """
  2. from __future__ import annotations
  3. from .add import Add
  4. from .mul import Mul, _keep_coeff
  5. from .power import Pow
  6. from .basic import Basic
  7. from .expr import Expr
  8. from .function import expand_power_exp
  9. from .sympify import sympify
  10. from .numbers import Rational, Integer, Number, I, equal_valued
  11. from .singleton import S
  12. from .sorting import default_sort_key, ordered
  13. from .symbol import Dummy
  14. from .traversal import preorder_traversal
  15. from .coreerrors import NonCommutativeExpression
  16. from .containers import Tuple, Dict
  17. from sympy.external.gmpy import SYMPY_INTS
  18. from sympy.utilities.iterables import (common_prefix, common_suffix,
  19. variations, iterable, is_sequence)
  20. from collections import defaultdict
  21. _eps = Dummy(positive=True)
  22. def _isnumber(i):
  23. return isinstance(i, (SYMPY_INTS, float)) or i.is_Number
  24. def _monotonic_sign(self):
  25. """Return the value closest to 0 that ``self`` may have if all symbols
  26. are signed and the result is uniformly the same sign for all values of symbols.
  27. If a symbol is only signed but not known to be an
  28. integer or the result is 0 then a symbol representative of the sign of self
  29. will be returned. Otherwise, None is returned if a) the sign could be positive
  30. or negative or b) self is not in one of the following forms:
  31. - L(x, y, ...) + A: a function linear in all symbols x, y, ... with an
  32. additive constant; if A is zero then the function can be a monomial whose
  33. sign is monotonic over the range of the variables, e.g. (x + 1)**3 if x is
  34. nonnegative.
  35. - A/L(x, y, ...) + B: the inverse of a function linear in all symbols x, y, ...
  36. that does not have a sign change from positive to negative for any set
  37. of values for the variables.
  38. - M(x, y, ...) + A: a monomial M whose factors are all signed and a constant, A.
  39. - A/M(x, y, ...) + B: the inverse of a monomial and constants A and B.
  40. - P(x): a univariate polynomial
  41. Examples
  42. ========
  43. >>> from sympy.core.exprtools import _monotonic_sign as F
  44. >>> from sympy import Dummy
  45. >>> nn = Dummy(integer=True, nonnegative=True)
  46. >>> p = Dummy(integer=True, positive=True)
  47. >>> p2 = Dummy(integer=True, positive=True)
  48. >>> F(nn + 1)
  49. 1
  50. >>> F(p - 1)
  51. _nneg
  52. >>> F(nn*p + 1)
  53. 1
  54. >>> F(p2*p + 1)
  55. 2
  56. >>> F(nn - 1) # could be negative, zero or positive
  57. """
  58. if not self.is_extended_real:
  59. return
  60. if (-self).is_Symbol:
  61. rv = _monotonic_sign(-self)
  62. return rv if rv is None else -rv
  63. if not self.is_Add and self.as_numer_denom()[1].is_number:
  64. s = self
  65. if s.is_prime:
  66. if s.is_odd:
  67. return Integer(3)
  68. else:
  69. return Integer(2)
  70. elif s.is_composite:
  71. if s.is_odd:
  72. return Integer(9)
  73. else:
  74. return Integer(4)
  75. elif s.is_positive:
  76. if s.is_even:
  77. if s.is_prime is False:
  78. return Integer(4)
  79. else:
  80. return Integer(2)
  81. elif s.is_integer:
  82. return S.One
  83. else:
  84. return _eps
  85. elif s.is_extended_negative:
  86. if s.is_even:
  87. return Integer(-2)
  88. elif s.is_integer:
  89. return S.NegativeOne
  90. else:
  91. return -_eps
  92. if s.is_zero or s.is_extended_nonpositive or s.is_extended_nonnegative:
  93. return S.Zero
  94. return None
  95. # univariate polynomial
  96. free = self.free_symbols
  97. if len(free) == 1:
  98. if self.is_polynomial():
  99. from sympy.polys.polytools import real_roots
  100. from sympy.polys.polyroots import roots
  101. from sympy.polys.polyerrors import PolynomialError
  102. x = free.pop()
  103. x0 = _monotonic_sign(x)
  104. if x0 in (_eps, -_eps):
  105. x0 = S.Zero
  106. if x0 is not None:
  107. d = self.diff(x)
  108. if d.is_number:
  109. currentroots = []
  110. else:
  111. try:
  112. currentroots = real_roots(d)
  113. except (PolynomialError, NotImplementedError):
  114. currentroots = [r for r in roots(d, x) if r.is_extended_real]
  115. y = self.subs(x, x0)
  116. if x.is_nonnegative and all(
  117. (r - x0).is_nonpositive for r in currentroots):
  118. if y.is_nonnegative and d.is_positive:
  119. if y:
  120. return y if y.is_positive else Dummy('pos', positive=True)
  121. else:
  122. return Dummy('nneg', nonnegative=True)
  123. if y.is_nonpositive and d.is_negative:
  124. if y:
  125. return y if y.is_negative else Dummy('neg', negative=True)
  126. else:
  127. return Dummy('npos', nonpositive=True)
  128. elif x.is_nonpositive and all(
  129. (r - x0).is_nonnegative for r in currentroots):
  130. if y.is_nonnegative and d.is_negative:
  131. if y:
  132. return Dummy('pos', positive=True)
  133. else:
  134. return Dummy('nneg', nonnegative=True)
  135. if y.is_nonpositive and d.is_positive:
  136. if y:
  137. return Dummy('neg', negative=True)
  138. else:
  139. return Dummy('npos', nonpositive=True)
  140. else:
  141. n, d = self.as_numer_denom()
  142. den = None
  143. if n.is_number:
  144. den = _monotonic_sign(d)
  145. elif not d.is_number:
  146. if _monotonic_sign(n) is not None:
  147. den = _monotonic_sign(d)
  148. if den is not None and (den.is_positive or den.is_negative):
  149. v = n*den
  150. if v.is_positive:
  151. return Dummy('pos', positive=True)
  152. elif v.is_nonnegative:
  153. return Dummy('nneg', nonnegative=True)
  154. elif v.is_negative:
  155. return Dummy('neg', negative=True)
  156. elif v.is_nonpositive:
  157. return Dummy('npos', nonpositive=True)
  158. return None
  159. # multivariate
  160. c, a = self.as_coeff_Add()
  161. v = None
  162. if not a.is_polynomial():
  163. # F/A or A/F where A is a number and F is a signed, rational monomial
  164. n, d = a.as_numer_denom()
  165. if not (n.is_number or d.is_number):
  166. return
  167. if (
  168. a.is_Mul or a.is_Pow) and \
  169. a.is_rational and \
  170. all(p.exp.is_Integer for p in a.atoms(Pow) if p.is_Pow) and \
  171. (a.is_positive or a.is_negative):
  172. v = S.One
  173. for ai in Mul.make_args(a):
  174. if ai.is_number:
  175. v *= ai
  176. continue
  177. reps = {}
  178. for x in ai.free_symbols:
  179. reps[x] = _monotonic_sign(x)
  180. if reps[x] is None:
  181. return
  182. v *= ai.subs(reps)
  183. elif c:
  184. # signed linear expression
  185. if not any(p for p in a.atoms(Pow) if not p.is_number) and (a.is_nonpositive or a.is_nonnegative):
  186. free = list(a.free_symbols)
  187. p = {}
  188. for i in free:
  189. v = _monotonic_sign(i)
  190. if v is None:
  191. return
  192. p[i] = v or (_eps if i.is_nonnegative else -_eps)
  193. v = a.xreplace(p)
  194. if v is not None:
  195. rv = v + c
  196. if v.is_nonnegative and rv.is_positive:
  197. return rv.subs(_eps, 0)
  198. if v.is_nonpositive and rv.is_negative:
  199. return rv.subs(_eps, 0)
  200. def decompose_power(expr: Expr) -> tuple[Expr, int]:
  201. """
  202. Decompose power into symbolic base and integer exponent.
  203. Examples
  204. ========
  205. >>> from sympy.core.exprtools import decompose_power
  206. >>> from sympy.abc import x, y
  207. >>> from sympy import exp
  208. >>> decompose_power(x)
  209. (x, 1)
  210. >>> decompose_power(x**2)
  211. (x, 2)
  212. >>> decompose_power(exp(2*y/3))
  213. (exp(y/3), 2)
  214. """
  215. base, exp = expr.as_base_exp()
  216. if exp.is_Number:
  217. if exp.is_Rational:
  218. if not exp.is_Integer:
  219. base = Pow(base, Rational(1, exp.q)) # type: ignore
  220. e = exp.p # type: ignore
  221. else:
  222. base, e = expr, 1
  223. else:
  224. exp, tail = exp.as_coeff_Mul(rational=True)
  225. if exp is S.NegativeOne:
  226. base, e = Pow(base, tail), -1
  227. elif exp is not S.One:
  228. # todo: after dropping python 3.7 support, use overload and Literal
  229. # in as_coeff_Mul to make exp Rational, and remove these 2 ignores
  230. tail = _keep_coeff(Rational(1, exp.q), tail) # type: ignore
  231. base, e = Pow(base, tail), exp.p # type: ignore
  232. else:
  233. base, e = expr, 1
  234. return base, e
  235. def decompose_power_rat(expr: Expr) -> tuple[Expr, Rational]:
  236. """
  237. Decompose power into symbolic base and rational exponent;
  238. if the exponent is not a Rational, then separate only the
  239. integer coefficient.
  240. Examples
  241. ========
  242. >>> from sympy.core.exprtools import decompose_power_rat
  243. >>> from sympy.abc import x
  244. >>> from sympy import sqrt, exp
  245. >>> decompose_power_rat(sqrt(x))
  246. (x, 1/2)
  247. >>> decompose_power_rat(exp(-3*x/2))
  248. (exp(x/2), -3)
  249. """
  250. base, exp = expr.as_base_exp()
  251. if not exp.is_Rational:
  252. base, exp_i = decompose_power(expr)
  253. exp = Integer(exp_i)
  254. return base, exp # type: ignore
  255. class Factors:
  256. """Efficient representation of ``f_1*f_2*...*f_n``."""
  257. __slots__ = ('factors', 'gens')
  258. def __init__(self, factors=None): # Factors
  259. """Initialize Factors from dict or expr.
  260. Examples
  261. ========
  262. >>> from sympy.core.exprtools import Factors
  263. >>> from sympy.abc import x
  264. >>> from sympy import I
  265. >>> e = 2*x**3
  266. >>> Factors(e)
  267. Factors({2: 1, x: 3})
  268. >>> Factors(e.as_powers_dict())
  269. Factors({2: 1, x: 3})
  270. >>> f = _
  271. >>> f.factors # underlying dictionary
  272. {2: 1, x: 3}
  273. >>> f.gens # base of each factor
  274. frozenset({2, x})
  275. >>> Factors(0)
  276. Factors({0: 1})
  277. >>> Factors(I)
  278. Factors({I: 1})
  279. Notes
  280. =====
  281. Although a dictionary can be passed, only minimal checking is
  282. performed: powers of -1 and I are made canonical.
  283. """
  284. if isinstance(factors, (SYMPY_INTS, float)):
  285. factors = S(factors)
  286. if isinstance(factors, Factors):
  287. factors = factors.factors.copy()
  288. elif factors in (None, S.One):
  289. factors = {}
  290. elif factors is S.Zero or factors == 0:
  291. factors = {S.Zero: S.One}
  292. elif isinstance(factors, Number):
  293. n = factors
  294. factors = {}
  295. if n < 0:
  296. factors[S.NegativeOne] = S.One
  297. n = -n
  298. if n is not S.One:
  299. if n.is_Float or n.is_Integer or n is S.Infinity:
  300. factors[n] = S.One
  301. elif n.is_Rational:
  302. # since we're processing Numbers, the denominator is
  303. # stored with a negative exponent; all other factors
  304. # are left .
  305. if n.p != 1:
  306. factors[Integer(n.p)] = S.One
  307. factors[Integer(n.q)] = S.NegativeOne
  308. else:
  309. raise ValueError('Expected Float|Rational|Integer, not %s' % n)
  310. elif isinstance(factors, Basic) and not factors.args:
  311. factors = {factors: S.One}
  312. elif isinstance(factors, Expr):
  313. c, nc = factors.args_cnc()
  314. i = c.count(I)
  315. for _ in range(i):
  316. c.remove(I)
  317. factors = dict(Mul._from_args(c).as_powers_dict())
  318. # Handle all rational Coefficients
  319. for f in list(factors.keys()):
  320. if isinstance(f, Rational) and not isinstance(f, Integer):
  321. p, q = Integer(f.p), Integer(f.q)
  322. factors[p] = (factors[p] if p in factors else S.Zero) + factors[f]
  323. factors[q] = (factors[q] if q in factors else S.Zero) - factors[f]
  324. factors.pop(f)
  325. if i:
  326. factors[I] = factors.get(I, S.Zero) + i
  327. if nc:
  328. factors[Mul(*nc, evaluate=False)] = S.One
  329. else:
  330. factors = factors.copy() # /!\ should be dict-like
  331. # tidy up -/+1 and I exponents if Rational
  332. handle = [k for k in factors if k is I or k in (-1, 1)]
  333. if handle:
  334. i1 = S.One
  335. for k in handle:
  336. if not _isnumber(factors[k]):
  337. continue
  338. i1 *= k**factors.pop(k)
  339. if i1 is not S.One:
  340. for a in i1.args if i1.is_Mul else [i1]: # at worst, -1.0*I*(-1)**e
  341. if a is S.NegativeOne:
  342. factors[a] = S.One
  343. elif a is I:
  344. factors[I] = S.One
  345. elif a.is_Pow:
  346. factors[a.base] = factors.get(a.base, S.Zero) + a.exp
  347. elif equal_valued(a, 1):
  348. factors[a] = S.One
  349. elif equal_valued(a, -1):
  350. factors[-a] = S.One
  351. factors[S.NegativeOne] = S.One
  352. else:
  353. raise ValueError('unexpected factor in i1: %s' % a)
  354. self.factors = factors
  355. keys = getattr(factors, 'keys', None)
  356. if keys is None:
  357. raise TypeError('expecting Expr or dictionary')
  358. self.gens = frozenset(keys())
  359. def __hash__(self): # Factors
  360. keys = tuple(ordered(self.factors.keys()))
  361. values = [self.factors[k] for k in keys]
  362. return hash((keys, values))
  363. def __repr__(self): # Factors
  364. return "Factors({%s})" % ', '.join(
  365. ['%s: %s' % (k, v) for k, v in ordered(self.factors.items())])
  366. @property
  367. def is_zero(self): # Factors
  368. """
  369. >>> from sympy.core.exprtools import Factors
  370. >>> Factors(0).is_zero
  371. True
  372. """
  373. f = self.factors
  374. return len(f) == 1 and S.Zero in f
  375. @property
  376. def is_one(self): # Factors
  377. """
  378. >>> from sympy.core.exprtools import Factors
  379. >>> Factors(1).is_one
  380. True
  381. """
  382. return not self.factors
  383. def as_expr(self): # Factors
  384. """Return the underlying expression.
  385. Examples
  386. ========
  387. >>> from sympy.core.exprtools import Factors
  388. >>> from sympy.abc import x, y
  389. >>> Factors((x*y**2).as_powers_dict()).as_expr()
  390. x*y**2
  391. """
  392. args = []
  393. for factor, exp in self.factors.items():
  394. if exp != 1:
  395. if isinstance(exp, Integer):
  396. b, e = factor.as_base_exp()
  397. e = _keep_coeff(exp, e)
  398. args.append(b**e)
  399. else:
  400. args.append(factor**exp)
  401. else:
  402. args.append(factor)
  403. return Mul(*args)
  404. def mul(self, other): # Factors
  405. """Return Factors of ``self * other``.
  406. Examples
  407. ========
  408. >>> from sympy.core.exprtools import Factors
  409. >>> from sympy.abc import x, y, z
  410. >>> a = Factors((x*y**2).as_powers_dict())
  411. >>> b = Factors((x*y/z).as_powers_dict())
  412. >>> a.mul(b)
  413. Factors({x: 2, y: 3, z: -1})
  414. >>> a*b
  415. Factors({x: 2, y: 3, z: -1})
  416. """
  417. if not isinstance(other, Factors):
  418. other = Factors(other)
  419. if any(f.is_zero for f in (self, other)):
  420. return Factors(S.Zero)
  421. factors = dict(self.factors)
  422. for factor, exp in other.factors.items():
  423. if factor in factors:
  424. exp = factors[factor] + exp
  425. if not exp:
  426. del factors[factor]
  427. continue
  428. factors[factor] = exp
  429. return Factors(factors)
  430. def normal(self, other):
  431. """Return ``self`` and ``other`` with ``gcd`` removed from each.
  432. The only differences between this and method ``div`` is that this
  433. is 1) optimized for the case when there are few factors in common and
  434. 2) this does not raise an error if ``other`` is zero.
  435. See Also
  436. ========
  437. div
  438. """
  439. if not isinstance(other, Factors):
  440. other = Factors(other)
  441. if other.is_zero:
  442. return (Factors(), Factors(S.Zero))
  443. if self.is_zero:
  444. return (Factors(S.Zero), Factors())
  445. self_factors = dict(self.factors)
  446. other_factors = dict(other.factors)
  447. for factor, self_exp in self.factors.items():
  448. try:
  449. other_exp = other.factors[factor]
  450. except KeyError:
  451. continue
  452. exp = self_exp - other_exp
  453. if not exp:
  454. del self_factors[factor]
  455. del other_factors[factor]
  456. elif _isnumber(exp):
  457. if exp > 0:
  458. self_factors[factor] = exp
  459. del other_factors[factor]
  460. else:
  461. del self_factors[factor]
  462. other_factors[factor] = -exp
  463. else:
  464. r = self_exp.extract_additively(other_exp)
  465. if r is not None:
  466. if r:
  467. self_factors[factor] = r
  468. del other_factors[factor]
  469. else: # should be handled already
  470. del self_factors[factor]
  471. del other_factors[factor]
  472. else:
  473. sc, sa = self_exp.as_coeff_Add()
  474. if sc:
  475. oc, oa = other_exp.as_coeff_Add()
  476. diff = sc - oc
  477. if diff > 0:
  478. self_factors[factor] -= oc
  479. other_exp = oa
  480. elif diff < 0:
  481. self_factors[factor] -= sc
  482. other_factors[factor] -= sc
  483. other_exp = oa - diff
  484. else:
  485. self_factors[factor] = sa
  486. other_exp = oa
  487. if other_exp:
  488. other_factors[factor] = other_exp
  489. else:
  490. del other_factors[factor]
  491. return Factors(self_factors), Factors(other_factors)
  492. def div(self, other): # Factors
  493. """Return ``self`` and ``other`` with ``gcd`` removed from each.
  494. This is optimized for the case when there are many factors in common.
  495. Examples
  496. ========
  497. >>> from sympy.core.exprtools import Factors
  498. >>> from sympy.abc import x, y, z
  499. >>> from sympy import S
  500. >>> a = Factors((x*y**2).as_powers_dict())
  501. >>> a.div(a)
  502. (Factors({}), Factors({}))
  503. >>> a.div(x*z)
  504. (Factors({y: 2}), Factors({z: 1}))
  505. The ``/`` operator only gives ``quo``:
  506. >>> a/x
  507. Factors({y: 2})
  508. Factors treats its factors as though they are all in the numerator, so
  509. if you violate this assumption the results will be correct but will
  510. not strictly correspond to the numerator and denominator of the ratio:
  511. >>> a.div(x/z)
  512. (Factors({y: 2}), Factors({z: -1}))
  513. Factors is also naive about bases: it does not attempt any denesting
  514. of Rational-base terms, for example the following does not become
  515. 2**(2*x)/2.
  516. >>> Factors(2**(2*x + 2)).div(S(8))
  517. (Factors({2: 2*x + 2}), Factors({8: 1}))
  518. factor_terms can clean up such Rational-bases powers:
  519. >>> from sympy import factor_terms
  520. >>> n, d = Factors(2**(2*x + 2)).div(S(8))
  521. >>> n.as_expr()/d.as_expr()
  522. 2**(2*x + 2)/8
  523. >>> factor_terms(_)
  524. 2**(2*x)/2
  525. """
  526. quo, rem = dict(self.factors), {}
  527. if not isinstance(other, Factors):
  528. other = Factors(other)
  529. if other.is_zero:
  530. raise ZeroDivisionError
  531. if self.is_zero:
  532. return (Factors(S.Zero), Factors())
  533. for factor, exp in other.factors.items():
  534. if factor in quo:
  535. d = quo[factor] - exp
  536. if _isnumber(d):
  537. if d <= 0:
  538. del quo[factor]
  539. if d >= 0:
  540. if d:
  541. quo[factor] = d
  542. continue
  543. exp = -d
  544. else:
  545. r = quo[factor].extract_additively(exp)
  546. if r is not None:
  547. if r:
  548. quo[factor] = r
  549. else: # should be handled already
  550. del quo[factor]
  551. else:
  552. other_exp = exp
  553. sc, sa = quo[factor].as_coeff_Add()
  554. if sc:
  555. oc, oa = other_exp.as_coeff_Add()
  556. diff = sc - oc
  557. if diff > 0:
  558. quo[factor] -= oc
  559. other_exp = oa
  560. elif diff < 0:
  561. quo[factor] -= sc
  562. other_exp = oa - diff
  563. else:
  564. quo[factor] = sa
  565. other_exp = oa
  566. if other_exp:
  567. rem[factor] = other_exp
  568. else:
  569. assert factor not in rem
  570. continue
  571. rem[factor] = exp
  572. return Factors(quo), Factors(rem)
  573. def quo(self, other): # Factors
  574. """Return numerator Factor of ``self / other``.
  575. Examples
  576. ========
  577. >>> from sympy.core.exprtools import Factors
  578. >>> from sympy.abc import x, y, z
  579. >>> a = Factors((x*y**2).as_powers_dict())
  580. >>> b = Factors((x*y/z).as_powers_dict())
  581. >>> a.quo(b) # same as a/b
  582. Factors({y: 1})
  583. """
  584. return self.div(other)[0]
  585. def rem(self, other): # Factors
  586. """Return denominator Factors of ``self / other``.
  587. Examples
  588. ========
  589. >>> from sympy.core.exprtools import Factors
  590. >>> from sympy.abc import x, y, z
  591. >>> a = Factors((x*y**2).as_powers_dict())
  592. >>> b = Factors((x*y/z).as_powers_dict())
  593. >>> a.rem(b)
  594. Factors({z: -1})
  595. >>> a.rem(a)
  596. Factors({})
  597. """
  598. return self.div(other)[1]
  599. def pow(self, other): # Factors
  600. """Return self raised to a non-negative integer power.
  601. Examples
  602. ========
  603. >>> from sympy.core.exprtools import Factors
  604. >>> from sympy.abc import x, y
  605. >>> a = Factors((x*y**2).as_powers_dict())
  606. >>> a**2
  607. Factors({x: 2, y: 4})
  608. """
  609. if isinstance(other, Factors):
  610. other = other.as_expr()
  611. if other.is_Integer:
  612. other = int(other)
  613. if isinstance(other, SYMPY_INTS) and other >= 0:
  614. factors = {}
  615. if other:
  616. for factor, exp in self.factors.items():
  617. factors[factor] = exp*other
  618. return Factors(factors)
  619. else:
  620. raise ValueError("expected non-negative integer, got %s" % other)
  621. def gcd(self, other): # Factors
  622. """Return Factors of ``gcd(self, other)``. The keys are
  623. the intersection of factors with the minimum exponent for
  624. each factor.
  625. Examples
  626. ========
  627. >>> from sympy.core.exprtools import Factors
  628. >>> from sympy.abc import x, y, z
  629. >>> a = Factors((x*y**2).as_powers_dict())
  630. >>> b = Factors((x*y/z).as_powers_dict())
  631. >>> a.gcd(b)
  632. Factors({x: 1, y: 1})
  633. """
  634. if not isinstance(other, Factors):
  635. other = Factors(other)
  636. if other.is_zero:
  637. return Factors(self.factors)
  638. factors = {}
  639. for factor, exp in self.factors.items():
  640. factor, exp = sympify(factor), sympify(exp)
  641. if factor in other.factors:
  642. lt = (exp - other.factors[factor]).is_negative
  643. if lt == True:
  644. factors[factor] = exp
  645. elif lt == False:
  646. factors[factor] = other.factors[factor]
  647. return Factors(factors)
  648. def lcm(self, other): # Factors
  649. """Return Factors of ``lcm(self, other)`` which are
  650. the union of factors with the maximum exponent for
  651. each factor.
  652. Examples
  653. ========
  654. >>> from sympy.core.exprtools import Factors
  655. >>> from sympy.abc import x, y, z
  656. >>> a = Factors((x*y**2).as_powers_dict())
  657. >>> b = Factors((x*y/z).as_powers_dict())
  658. >>> a.lcm(b)
  659. Factors({x: 1, y: 2, z: -1})
  660. """
  661. if not isinstance(other, Factors):
  662. other = Factors(other)
  663. if any(f.is_zero for f in (self, other)):
  664. return Factors(S.Zero)
  665. factors = dict(self.factors)
  666. for factor, exp in other.factors.items():
  667. if factor in factors:
  668. exp = max(exp, factors[factor])
  669. factors[factor] = exp
  670. return Factors(factors)
  671. def __mul__(self, other): # Factors
  672. return self.mul(other)
  673. def __divmod__(self, other): # Factors
  674. return self.div(other)
  675. def __truediv__(self, other): # Factors
  676. return self.quo(other)
  677. def __mod__(self, other): # Factors
  678. return self.rem(other)
  679. def __pow__(self, other): # Factors
  680. return self.pow(other)
  681. def __eq__(self, other): # Factors
  682. if not isinstance(other, Factors):
  683. other = Factors(other)
  684. return self.factors == other.factors
  685. def __ne__(self, other): # Factors
  686. return not self == other
  687. class Term:
  688. """Efficient representation of ``coeff*(numer/denom)``. """
  689. __slots__ = ('coeff', 'numer', 'denom')
  690. def __init__(self, term, numer=None, denom=None): # Term
  691. if numer is None and denom is None:
  692. if not term.is_commutative:
  693. raise NonCommutativeExpression(
  694. 'commutative expression expected')
  695. coeff, factors = term.as_coeff_mul()
  696. numer, denom = defaultdict(int), defaultdict(int)
  697. for factor in factors:
  698. base, exp = decompose_power(factor)
  699. if base.is_Add:
  700. cont, base = base.primitive()
  701. coeff *= cont**exp
  702. if exp > 0:
  703. numer[base] += exp
  704. else:
  705. denom[base] += -exp
  706. numer = Factors(numer)
  707. denom = Factors(denom)
  708. else:
  709. coeff = term
  710. if numer is None:
  711. numer = Factors()
  712. if denom is None:
  713. denom = Factors()
  714. self.coeff = coeff
  715. self.numer = numer
  716. self.denom = denom
  717. def __hash__(self): # Term
  718. return hash((self.coeff, self.numer, self.denom))
  719. def __repr__(self): # Term
  720. return "Term(%s, %s, %s)" % (self.coeff, self.numer, self.denom)
  721. def as_expr(self): # Term
  722. return self.coeff*(self.numer.as_expr()/self.denom.as_expr())
  723. def mul(self, other): # Term
  724. coeff = self.coeff*other.coeff
  725. numer = self.numer.mul(other.numer)
  726. denom = self.denom.mul(other.denom)
  727. numer, denom = numer.normal(denom)
  728. return Term(coeff, numer, denom)
  729. def inv(self): # Term
  730. return Term(1/self.coeff, self.denom, self.numer)
  731. def quo(self, other): # Term
  732. return self.mul(other.inv())
  733. def pow(self, other): # Term
  734. if other < 0:
  735. return self.inv().pow(-other)
  736. else:
  737. return Term(self.coeff ** other,
  738. self.numer.pow(other),
  739. self.denom.pow(other))
  740. def gcd(self, other): # Term
  741. return Term(self.coeff.gcd(other.coeff),
  742. self.numer.gcd(other.numer),
  743. self.denom.gcd(other.denom))
  744. def lcm(self, other): # Term
  745. return Term(self.coeff.lcm(other.coeff),
  746. self.numer.lcm(other.numer),
  747. self.denom.lcm(other.denom))
  748. def __mul__(self, other): # Term
  749. if isinstance(other, Term):
  750. return self.mul(other)
  751. else:
  752. return NotImplemented
  753. def __truediv__(self, other): # Term
  754. if isinstance(other, Term):
  755. return self.quo(other)
  756. else:
  757. return NotImplemented
  758. def __pow__(self, other): # Term
  759. if isinstance(other, SYMPY_INTS):
  760. return self.pow(other)
  761. else:
  762. return NotImplemented
  763. def __eq__(self, other): # Term
  764. return (self.coeff == other.coeff and
  765. self.numer == other.numer and
  766. self.denom == other.denom)
  767. def __ne__(self, other): # Term
  768. return not self == other
  769. def _gcd_terms(terms, isprimitive=False, fraction=True):
  770. """Helper function for :func:`gcd_terms`.
  771. Parameters
  772. ==========
  773. isprimitive : boolean, optional
  774. If ``isprimitive`` is True then the call to primitive
  775. for an Add will be skipped. This is useful when the
  776. content has already been extracted.
  777. fraction : boolean, optional
  778. If ``fraction`` is True then the expression will appear over a common
  779. denominator, the lcm of all term denominators.
  780. """
  781. if isinstance(terms, Basic) and not isinstance(terms, Tuple):
  782. terms = Add.make_args(terms)
  783. terms = list(map(Term, [t for t in terms if t]))
  784. # there is some simplification that may happen if we leave this
  785. # here rather than duplicate it before the mapping of Term onto
  786. # the terms
  787. if len(terms) == 0:
  788. return S.Zero, S.Zero, S.One
  789. if len(terms) == 1:
  790. cont = terms[0].coeff
  791. numer = terms[0].numer.as_expr()
  792. denom = terms[0].denom.as_expr()
  793. else:
  794. cont = terms[0]
  795. for term in terms[1:]:
  796. cont = cont.gcd(term)
  797. for i, term in enumerate(terms):
  798. terms[i] = term.quo(cont)
  799. if fraction:
  800. denom = terms[0].denom
  801. for term in terms[1:]:
  802. denom = denom.lcm(term.denom)
  803. numers = []
  804. for term in terms:
  805. numer = term.numer.mul(denom.quo(term.denom))
  806. numers.append(term.coeff*numer.as_expr())
  807. else:
  808. numers = [t.as_expr() for t in terms]
  809. denom = Term(S.One).numer
  810. cont = cont.as_expr()
  811. numer = Add(*numers)
  812. denom = denom.as_expr()
  813. if not isprimitive and numer.is_Add:
  814. _cont, numer = numer.primitive()
  815. cont *= _cont
  816. return cont, numer, denom
  817. def gcd_terms(terms, isprimitive=False, clear=True, fraction=True):
  818. """Compute the GCD of ``terms`` and put them together.
  819. Parameters
  820. ==========
  821. terms : Expr
  822. Can be an expression or a non-Basic sequence of expressions
  823. which will be handled as though they are terms from a sum.
  824. isprimitive : bool, optional
  825. If ``isprimitive`` is True the _gcd_terms will not run the primitive
  826. method on the terms.
  827. clear : bool, optional
  828. It controls the removal of integers from the denominator of an Add
  829. expression. When True (default), all numerical denominator will be cleared;
  830. when False the denominators will be cleared only if all terms had numerical
  831. denominators other than 1.
  832. fraction : bool, optional
  833. When True (default), will put the expression over a common
  834. denominator.
  835. Examples
  836. ========
  837. >>> from sympy import gcd_terms
  838. >>> from sympy.abc import x, y
  839. >>> gcd_terms((x + 1)**2*y + (x + 1)*y**2)
  840. y*(x + 1)*(x + y + 1)
  841. >>> gcd_terms(x/2 + 1)
  842. (x + 2)/2
  843. >>> gcd_terms(x/2 + 1, clear=False)
  844. x/2 + 1
  845. >>> gcd_terms(x/2 + y/2, clear=False)
  846. (x + y)/2
  847. >>> gcd_terms(x/2 + 1/x)
  848. (x**2 + 2)/(2*x)
  849. >>> gcd_terms(x/2 + 1/x, fraction=False)
  850. (x + 2/x)/2
  851. >>> gcd_terms(x/2 + 1/x, fraction=False, clear=False)
  852. x/2 + 1/x
  853. >>> gcd_terms(x/2/y + 1/x/y)
  854. (x**2 + 2)/(2*x*y)
  855. >>> gcd_terms(x/2/y + 1/x/y, clear=False)
  856. (x**2/2 + 1)/(x*y)
  857. >>> gcd_terms(x/2/y + 1/x/y, clear=False, fraction=False)
  858. (x/2 + 1/x)/y
  859. The ``clear`` flag was ignored in this case because the returned
  860. expression was a rational expression, not a simple sum.
  861. See Also
  862. ========
  863. factor_terms, sympy.polys.polytools.terms_gcd
  864. """
  865. def mask(terms):
  866. """replace nc portions of each term with a unique Dummy symbols
  867. and return the replacements to restore them"""
  868. args = [(a, []) if a.is_commutative else a.args_cnc() for a in terms]
  869. reps = []
  870. for i, (c, nc) in enumerate(args):
  871. if nc:
  872. nc = Mul(*nc)
  873. d = Dummy()
  874. reps.append((d, nc))
  875. c.append(d)
  876. args[i] = Mul(*c)
  877. else:
  878. args[i] = c
  879. return args, dict(reps)
  880. isadd = isinstance(terms, Add)
  881. addlike = isadd or not isinstance(terms, Basic) and \
  882. is_sequence(terms, include=set) and \
  883. not isinstance(terms, Dict)
  884. if addlike:
  885. if isadd: # i.e. an Add
  886. terms = list(terms.args)
  887. else:
  888. terms = sympify(terms)
  889. terms, reps = mask(terms)
  890. cont, numer, denom = _gcd_terms(terms, isprimitive, fraction)
  891. numer = numer.xreplace(reps)
  892. coeff, factors = cont.as_coeff_Mul()
  893. if not clear:
  894. c, _coeff = coeff.as_coeff_Mul()
  895. if not c.is_Integer and not clear and numer.is_Add:
  896. n, d = c.as_numer_denom()
  897. _numer = numer/d
  898. if any(a.as_coeff_Mul()[0].is_Integer
  899. for a in _numer.args):
  900. numer = _numer
  901. coeff = n*_coeff
  902. return _keep_coeff(coeff, factors*numer/denom, clear=clear)
  903. if not isinstance(terms, Basic):
  904. return terms
  905. if terms.is_Atom:
  906. return terms
  907. if terms.is_Mul:
  908. c, args = terms.as_coeff_mul()
  909. return _keep_coeff(c, Mul(*[gcd_terms(i, isprimitive, clear, fraction)
  910. for i in args]), clear=clear)
  911. def handle(a):
  912. # don't treat internal args like terms of an Add
  913. if not isinstance(a, Expr):
  914. if isinstance(a, Basic):
  915. if not a.args:
  916. return a
  917. return a.func(*[handle(i) for i in a.args])
  918. return type(a)([handle(i) for i in a])
  919. return gcd_terms(a, isprimitive, clear, fraction)
  920. if isinstance(terms, Dict):
  921. return Dict(*[(k, handle(v)) for k, v in terms.args])
  922. return terms.func(*[handle(i) for i in terms.args])
  923. def _factor_sum_int(expr, **kwargs):
  924. """Return Sum or Integral object with factors that are not
  925. in the wrt variables removed. In cases where there are additive
  926. terms in the function of the object that are independent, the
  927. object will be separated into two objects.
  928. Examples
  929. ========
  930. >>> from sympy import Sum, factor_terms
  931. >>> from sympy.abc import x, y
  932. >>> factor_terms(Sum(x + y, (x, 1, 3)))
  933. y*Sum(1, (x, 1, 3)) + Sum(x, (x, 1, 3))
  934. >>> factor_terms(Sum(x*y, (x, 1, 3)))
  935. y*Sum(x, (x, 1, 3))
  936. Notes
  937. =====
  938. If a function in the summand or integrand is replaced
  939. with a symbol, then this simplification should not be
  940. done or else an incorrect result will be obtained when
  941. the symbol is replaced with an expression that depends
  942. on the variables of summation/integration:
  943. >>> eq = Sum(y, (x, 1, 3))
  944. >>> factor_terms(eq).subs(y, x).doit()
  945. 3*x
  946. >>> eq.subs(y, x).doit()
  947. 6
  948. """
  949. result = expr.function
  950. if result == 0:
  951. return S.Zero
  952. limits = expr.limits
  953. # get the wrt variables
  954. wrt = {i.args[0] for i in limits}
  955. # factor out any common terms that are independent of wrt
  956. f = factor_terms(result, **kwargs)
  957. i, d = f.as_independent(*wrt)
  958. if isinstance(f, Add):
  959. return i * expr.func(1, *limits) + expr.func(d, *limits)
  960. else:
  961. return i * expr.func(d, *limits)
  962. def factor_terms(expr: Expr | complex, radical=False, clear=False, fraction=False, sign=True) -> Expr:
  963. """Remove common factors from terms in all arguments without
  964. changing the underlying structure of the expr. No expansion or
  965. simplification (and no processing of non-commutatives) is performed.
  966. Parameters
  967. ==========
  968. radical: bool, optional
  969. If radical=True then a radical common to all terms will be factored
  970. out of any Add sub-expressions of the expr.
  971. clear : bool, optional
  972. If clear=False (default) then coefficients will not be separated
  973. from a single Add if they can be distributed to leave one or more
  974. terms with integer coefficients.
  975. fraction : bool, optional
  976. If fraction=True (default is False) then a common denominator will be
  977. constructed for the expression.
  978. sign : bool, optional
  979. If sign=True (default) then even if the only factor in common is a -1,
  980. it will be factored out of the expression.
  981. Examples
  982. ========
  983. >>> from sympy import factor_terms, Symbol
  984. >>> from sympy.abc import x, y
  985. >>> factor_terms(x + x*(2 + 4*y)**3)
  986. x*(8*(2*y + 1)**3 + 1)
  987. >>> A = Symbol('A', commutative=False)
  988. >>> factor_terms(x*A + x*A + x*y*A)
  989. x*(y*A + 2*A)
  990. When ``clear`` is False, a rational will only be factored out of an
  991. Add expression if all terms of the Add have coefficients that are
  992. fractions:
  993. >>> factor_terms(x/2 + 1, clear=False)
  994. x/2 + 1
  995. >>> factor_terms(x/2 + 1, clear=True)
  996. (x + 2)/2
  997. If a -1 is all that can be factored out, to *not* factor it out, the
  998. flag ``sign`` must be False:
  999. >>> factor_terms(-x - y)
  1000. -(x + y)
  1001. >>> factor_terms(-x - y, sign=False)
  1002. -x - y
  1003. >>> factor_terms(-2*x - 2*y, sign=False)
  1004. -2*(x + y)
  1005. See Also
  1006. ========
  1007. gcd_terms, sympy.polys.polytools.terms_gcd
  1008. """
  1009. def do(expr):
  1010. from sympy.concrete.summations import Sum
  1011. from sympy.integrals.integrals import Integral
  1012. is_iterable = iterable(expr)
  1013. if not isinstance(expr, Basic) or expr.is_Atom:
  1014. if is_iterable:
  1015. return type(expr)([do(i) for i in expr])
  1016. return expr
  1017. if expr.is_Pow or expr.is_Function or \
  1018. is_iterable or not hasattr(expr, 'args_cnc'):
  1019. args = expr.args
  1020. newargs = tuple([do(i) for i in args])
  1021. if newargs == args:
  1022. return expr
  1023. return expr.func(*newargs)
  1024. if isinstance(expr, (Sum, Integral)):
  1025. return _factor_sum_int(expr,
  1026. radical=radical, clear=clear,
  1027. fraction=fraction, sign=sign)
  1028. cont, p = expr.as_content_primitive(radical=radical, clear=clear)
  1029. if p.is_Add:
  1030. list_args = [do(a) for a in Add.make_args(p)]
  1031. # get a common negative (if there) which gcd_terms does not remove
  1032. if not any(a.as_coeff_Mul()[0].extract_multiplicatively(-1) is None
  1033. for a in list_args):
  1034. cont = -cont
  1035. list_args = [-a for a in list_args]
  1036. # watch out for exp(-(x+2)) which gcd_terms will change to exp(-x-2)
  1037. special = {}
  1038. for i, a in enumerate(list_args):
  1039. b, e = a.as_base_exp()
  1040. if e.is_Mul and e != Mul(*e.args):
  1041. list_args[i] = Dummy()
  1042. special[list_args[i]] = a
  1043. # rebuild p not worrying about the order which gcd_terms will fix
  1044. p = Add._from_args(list_args)
  1045. p = gcd_terms(p,
  1046. isprimitive=True,
  1047. clear=clear,
  1048. fraction=fraction).xreplace(special)
  1049. elif p.args:
  1050. p = p.func(
  1051. *[do(a) for a in p.args])
  1052. rv = _keep_coeff(cont, p, clear=clear, sign=sign)
  1053. return rv
  1054. expr2 = sympify(expr)
  1055. return do(expr2)
  1056. def _mask_nc(eq, name=None):
  1057. """
  1058. Return ``eq`` with non-commutative objects replaced with Dummy
  1059. symbols. A dictionary that can be used to restore the original
  1060. values is returned: if it is None, the expression is noncommutative
  1061. and cannot be made commutative. The third value returned is a list
  1062. of any non-commutative symbols that appear in the returned equation.
  1063. Explanation
  1064. ===========
  1065. All non-commutative objects other than Symbols are replaced with
  1066. a non-commutative Symbol. Identical objects will be identified
  1067. by identical symbols.
  1068. If there is only 1 non-commutative object in an expression it will
  1069. be replaced with a commutative symbol. Otherwise, the non-commutative
  1070. entities are retained and the calling routine should handle
  1071. replacements in this case since some care must be taken to keep
  1072. track of the ordering of symbols when they occur within Muls.
  1073. Parameters
  1074. ==========
  1075. name : str
  1076. ``name``, if given, is the name that will be used with numbered Dummy
  1077. variables that will replace the non-commutative objects and is mainly
  1078. used for doctesting purposes.
  1079. Examples
  1080. ========
  1081. >>> from sympy.physics.secondquant import Commutator, NO, F, Fd
  1082. >>> from sympy import symbols
  1083. >>> from sympy.core.exprtools import _mask_nc
  1084. >>> from sympy.abc import x, y
  1085. >>> A, B, C = symbols('A,B,C', commutative=False)
  1086. One nc-symbol:
  1087. >>> _mask_nc(A**2 - x**2, 'd')
  1088. (_d0**2 - x**2, {_d0: A}, [])
  1089. Multiple nc-symbols:
  1090. >>> _mask_nc(A**2 - B**2, 'd')
  1091. (A**2 - B**2, {}, [A, B])
  1092. An nc-object with nc-symbols but no others outside of it:
  1093. >>> _mask_nc(1 + x*Commutator(A, B), 'd')
  1094. (_d0*x + 1, {_d0: Commutator(A, B)}, [])
  1095. >>> _mask_nc(NO(Fd(x)*F(y)), 'd')
  1096. (_d0, {_d0: NO(CreateFermion(x)*AnnihilateFermion(y))}, [])
  1097. Multiple nc-objects:
  1098. >>> eq = x*Commutator(A, B) + x*Commutator(A, C)*Commutator(A, B)
  1099. >>> _mask_nc(eq, 'd')
  1100. (x*_d0 + x*_d1*_d0, {_d0: Commutator(A, B), _d1: Commutator(A, C)}, [_d0, _d1])
  1101. Multiple nc-objects and nc-symbols:
  1102. >>> eq = A*Commutator(A, B) + B*Commutator(A, C)
  1103. >>> _mask_nc(eq, 'd')
  1104. (A*_d0 + B*_d1, {_d0: Commutator(A, B), _d1: Commutator(A, C)}, [_d0, _d1, A, B])
  1105. """
  1106. name = name or 'mask'
  1107. # Make Dummy() append sequential numbers to the name
  1108. def numbered_names():
  1109. i = 0
  1110. while True:
  1111. yield name + str(i)
  1112. i += 1
  1113. names = numbered_names()
  1114. def Dummy(*args, **kwargs):
  1115. from .symbol import Dummy
  1116. return Dummy(next(names), *args, **kwargs)
  1117. expr = eq
  1118. if expr.is_commutative:
  1119. return eq, {}, []
  1120. # identify nc-objects; symbols and other
  1121. rep = []
  1122. nc_obj = set()
  1123. nc_syms = set()
  1124. pot = preorder_traversal(expr, keys=default_sort_key)
  1125. for a in pot:
  1126. if any(a == r[0] for r in rep):
  1127. pot.skip()
  1128. elif not a.is_commutative:
  1129. if a.is_symbol:
  1130. nc_syms.add(a)
  1131. pot.skip()
  1132. elif not (a.is_Add or a.is_Mul or a.is_Pow):
  1133. nc_obj.add(a)
  1134. pot.skip()
  1135. # If there is only one nc symbol or object, it can be factored regularly
  1136. # but polys is going to complain, so replace it with a Dummy.
  1137. if len(nc_obj) == 1 and not nc_syms:
  1138. rep.append((nc_obj.pop(), Dummy()))
  1139. elif len(nc_syms) == 1 and not nc_obj:
  1140. rep.append((nc_syms.pop(), Dummy()))
  1141. # Any remaining nc-objects will be replaced with an nc-Dummy and
  1142. # identified as an nc-Symbol to watch out for
  1143. nc_obj = sorted(nc_obj, key=default_sort_key)
  1144. for n in nc_obj:
  1145. nc = Dummy(commutative=False)
  1146. rep.append((n, nc))
  1147. nc_syms.add(nc)
  1148. expr = expr.subs(rep)
  1149. nc_syms = list(nc_syms)
  1150. nc_syms.sort(key=default_sort_key)
  1151. return expr, {v: k for k, v in rep}, nc_syms
  1152. def factor_nc(expr):
  1153. """Return the factored form of ``expr`` while handling non-commutative
  1154. expressions.
  1155. Examples
  1156. ========
  1157. >>> from sympy import factor_nc, Symbol
  1158. >>> from sympy.abc import x
  1159. >>> A = Symbol('A', commutative=False)
  1160. >>> B = Symbol('B', commutative=False)
  1161. >>> factor_nc((x**2 + 2*A*x + A**2).expand())
  1162. (x + A)**2
  1163. >>> factor_nc(((x + A)*(x + B)).expand())
  1164. (x + A)*(x + B)
  1165. """
  1166. expr = sympify(expr)
  1167. if not isinstance(expr, Expr) or not expr.args:
  1168. return expr
  1169. if not expr.is_Add:
  1170. return expr.func(*[factor_nc(a) for a in expr.args])
  1171. expr = expr.func(*[expand_power_exp(i) for i in expr.args])
  1172. from sympy.polys.polytools import gcd, factor
  1173. expr, rep, nc_symbols = _mask_nc(expr)
  1174. if rep:
  1175. return factor(expr).subs(rep)
  1176. else:
  1177. args = [a.args_cnc() for a in Add.make_args(expr)]
  1178. c = g = l = r = S.One
  1179. hit = False
  1180. # find any commutative gcd term
  1181. for i, a in enumerate(args):
  1182. if i == 0:
  1183. c = Mul._from_args(a[0])
  1184. elif a[0]:
  1185. c = gcd(c, Mul._from_args(a[0]))
  1186. else:
  1187. c = S.One
  1188. if c is not S.One:
  1189. hit = True
  1190. c, g = c.as_coeff_Mul()
  1191. if g is not S.One:
  1192. for i, (cc, _) in enumerate(args):
  1193. cc = list(Mul.make_args(Mul._from_args(list(cc))/g))
  1194. args[i][0] = cc
  1195. for i, (cc, _) in enumerate(args):
  1196. if cc:
  1197. cc[0] = cc[0]/c
  1198. else:
  1199. cc = [1/c]
  1200. args[i][0] = cc
  1201. # find any noncommutative common prefix
  1202. for i, a in enumerate(args):
  1203. if i == 0:
  1204. n = a[1][:]
  1205. else:
  1206. n = common_prefix(n, a[1])
  1207. if not n:
  1208. # is there a power that can be extracted?
  1209. if not args[0][1]:
  1210. break
  1211. b, e = args[0][1][0].as_base_exp()
  1212. ok = False
  1213. if e.is_Integer:
  1214. for t in args:
  1215. if not t[1]:
  1216. break
  1217. bt, et = t[1][0].as_base_exp()
  1218. if et.is_Integer and bt == b:
  1219. e = min(e, et)
  1220. else:
  1221. break
  1222. else:
  1223. ok = hit = True
  1224. l = b**e
  1225. il = b**-e
  1226. for _ in args:
  1227. _[1][0] = il*_[1][0]
  1228. break
  1229. if not ok:
  1230. break
  1231. else:
  1232. hit = True
  1233. lenn = len(n)
  1234. l = Mul(*n)
  1235. for _ in args:
  1236. _[1] = _[1][lenn:]
  1237. # find any noncommutative common suffix
  1238. for i, a in enumerate(args):
  1239. if i == 0:
  1240. n = a[1][:]
  1241. else:
  1242. n = common_suffix(n, a[1])
  1243. if not n:
  1244. # is there a power that can be extracted?
  1245. if not args[0][1]:
  1246. break
  1247. b, e = args[0][1][-1].as_base_exp()
  1248. ok = False
  1249. if e.is_Integer:
  1250. for t in args:
  1251. if not t[1]:
  1252. break
  1253. bt, et = t[1][-1].as_base_exp()
  1254. if et.is_Integer and bt == b:
  1255. e = min(e, et)
  1256. else:
  1257. break
  1258. else:
  1259. ok = hit = True
  1260. r = b**e
  1261. il = b**-e
  1262. for _ in args:
  1263. _[1][-1] = _[1][-1]*il
  1264. break
  1265. if not ok:
  1266. break
  1267. else:
  1268. hit = True
  1269. lenn = len(n)
  1270. r = Mul(*n)
  1271. for _ in args:
  1272. _[1] = _[1][:len(_[1]) - lenn]
  1273. if hit:
  1274. mid = Add(*[Mul(*cc)*Mul(*nc) for cc, nc in args])
  1275. else:
  1276. mid = expr
  1277. from sympy.simplify.powsimp import powsimp
  1278. # sort the symbols so the Dummys would appear in the same
  1279. # order as the original symbols, otherwise you may introduce
  1280. # a factor of -1, e.g. A**2 - B**2) -- {A:y, B:x} --> y**2 - x**2
  1281. # and the former factors into two terms, (A - B)*(A + B) while the
  1282. # latter factors into 3 terms, (-1)*(x - y)*(x + y)
  1283. rep1 = [(n, Dummy()) for n in sorted(nc_symbols, key=default_sort_key)]
  1284. unrep1 = [(v, k) for k, v in rep1]
  1285. unrep1.reverse()
  1286. new_mid, r2, _ = _mask_nc(mid.subs(rep1))
  1287. new_mid = powsimp(factor(new_mid))
  1288. new_mid = new_mid.subs(r2).subs(unrep1)
  1289. if new_mid.is_Pow:
  1290. return _keep_coeff(c, g*l*new_mid*r)
  1291. if new_mid.is_Mul:
  1292. def _pemexpand(expr):
  1293. "Expand with the minimal set of hints necessary to check the result."
  1294. return expr.expand(deep=True, mul=True, power_exp=True,
  1295. power_base=False, basic=False, multinomial=True, log=False)
  1296. # XXX TODO there should be a way to inspect what order the terms
  1297. # must be in and just select the plausible ordering without
  1298. # checking permutations
  1299. cfac = []
  1300. ncfac = []
  1301. for f in new_mid.args:
  1302. if f.is_commutative:
  1303. cfac.append(f)
  1304. else:
  1305. b, e = f.as_base_exp()
  1306. if e.is_Integer:
  1307. ncfac.extend([b]*e)
  1308. else:
  1309. ncfac.append(f)
  1310. pre_mid = g*Mul(*cfac)*l
  1311. target = _pemexpand(expr/c)
  1312. for s in variations(ncfac, len(ncfac)):
  1313. ok = pre_mid*Mul(*s)*r
  1314. if _pemexpand(ok) == target:
  1315. return _keep_coeff(c, ok)
  1316. # mid was an Add that didn't factor successfully
  1317. return _keep_coeff(c, g*l*mid*r)