simplify.py 72 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164
  1. from __future__ import annotations
  2. from typing import overload
  3. from collections import defaultdict
  4. from sympy.concrete.products import Product
  5. from sympy.concrete.summations import Sum
  6. from sympy.core import (Basic, S, Add, Mul, Pow, Symbol, sympify,
  7. expand_func, Function, Dummy, Expr, factor_terms,
  8. expand_power_exp, Eq)
  9. from sympy.core.exprtools import factor_nc
  10. from sympy.core.parameters import global_parameters
  11. from sympy.core.function import (expand_log, count_ops, _mexpand,
  12. nfloat, expand_mul, expand)
  13. from sympy.core.numbers import Float, I, pi, Rational, equal_valued
  14. from sympy.core.relational import Relational
  15. from sympy.core.rules import Transform
  16. from sympy.core.sorting import ordered
  17. from sympy.core.sympify import _sympify
  18. from sympy.core.traversal import bottom_up as _bottom_up, walk as _walk
  19. from sympy.functions import gamma, exp, sqrt, log, exp_polar, re
  20. from sympy.functions.combinatorial.factorials import CombinatorialFunction
  21. from sympy.functions.elementary.complexes import unpolarify, Abs, sign
  22. from sympy.functions.elementary.exponential import ExpBase
  23. from sympy.functions.elementary.hyperbolic import HyperbolicFunction
  24. from sympy.functions.elementary.integers import ceiling
  25. from sympy.functions.elementary.piecewise import (Piecewise, piecewise_fold,
  26. piecewise_simplify)
  27. from sympy.functions.elementary.trigonometric import TrigonometricFunction
  28. from sympy.functions.special.bessel import (BesselBase, besselj, besseli,
  29. besselk, bessely, jn)
  30. from sympy.functions.special.tensor_functions import KroneckerDelta
  31. from sympy.integrals.integrals import Integral
  32. from sympy.logic.boolalg import Boolean
  33. from sympy.matrices.expressions import (MatrixExpr, MatAdd, MatMul,
  34. MatPow, MatrixSymbol)
  35. from sympy.polys import together, cancel, factor
  36. from sympy.polys.numberfields.minpoly import _is_sum_surds, _minimal_polynomial_sq
  37. from sympy.sets.sets import Set
  38. from sympy.simplify.combsimp import combsimp
  39. from sympy.simplify.cse_opts import sub_pre, sub_post
  40. from sympy.simplify.hyperexpand import hyperexpand
  41. from sympy.simplify.powsimp import powsimp
  42. from sympy.simplify.radsimp import radsimp, fraction, collect_abs
  43. from sympy.simplify.sqrtdenest import sqrtdenest
  44. from sympy.simplify.trigsimp import trigsimp, exptrigsimp
  45. from sympy.utilities.decorator import deprecated
  46. from sympy.utilities.iterables import has_variety, sift, subsets, iterable
  47. from sympy.utilities.misc import as_int
  48. import mpmath
  49. def separatevars(expr, symbols=[], dict=False, force=False):
  50. """
  51. Separates variables in an expression, if possible. By
  52. default, it separates with respect to all symbols in an
  53. expression and collects constant coefficients that are
  54. independent of symbols.
  55. Explanation
  56. ===========
  57. If ``dict=True`` then the separated terms will be returned
  58. in a dictionary keyed to their corresponding symbols.
  59. By default, all symbols in the expression will appear as
  60. keys; if symbols are provided, then all those symbols will
  61. be used as keys, and any terms in the expression containing
  62. other symbols or non-symbols will be returned keyed to the
  63. string 'coeff'. (Passing None for symbols will return the
  64. expression in a dictionary keyed to 'coeff'.)
  65. If ``force=True``, then bases of powers will be separated regardless
  66. of assumptions on the symbols involved.
  67. Notes
  68. =====
  69. The order of the factors is determined by Mul, so that the
  70. separated expressions may not necessarily be grouped together.
  71. Although factoring is necessary to separate variables in some
  72. expressions, it is not necessary in all cases, so one should not
  73. count on the returned factors being factored.
  74. Examples
  75. ========
  76. >>> from sympy.abc import x, y, z, alpha
  77. >>> from sympy import separatevars, sin
  78. >>> separatevars((x*y)**y)
  79. (x*y)**y
  80. >>> separatevars((x*y)**y, force=True)
  81. x**y*y**y
  82. >>> e = 2*x**2*z*sin(y)+2*z*x**2
  83. >>> separatevars(e)
  84. 2*x**2*z*(sin(y) + 1)
  85. >>> separatevars(e, symbols=(x, y), dict=True)
  86. {'coeff': 2*z, x: x**2, y: sin(y) + 1}
  87. >>> separatevars(e, [x, y, alpha], dict=True)
  88. {'coeff': 2*z, alpha: 1, x: x**2, y: sin(y) + 1}
  89. If the expression is not really separable, or is only partially
  90. separable, separatevars will do the best it can to separate it
  91. by using factoring.
  92. >>> separatevars(x + x*y - 3*x**2)
  93. -x*(3*x - y - 1)
  94. If the expression is not separable then expr is returned unchanged
  95. or (if dict=True) then None is returned.
  96. >>> eq = 2*x + y*sin(x)
  97. >>> separatevars(eq) == eq
  98. True
  99. >>> separatevars(2*x + y*sin(x), symbols=(x, y), dict=True) is None
  100. True
  101. """
  102. expr = sympify(expr)
  103. if dict:
  104. return _separatevars_dict(_separatevars(expr, force), symbols)
  105. else:
  106. return _separatevars(expr, force)
  107. def _separatevars(expr, force):
  108. if isinstance(expr, Abs):
  109. arg = expr.args[0]
  110. if arg.is_Mul and not arg.is_number:
  111. s = separatevars(arg, dict=True, force=force)
  112. if s is not None:
  113. return Mul(*map(expr.func, s.values()))
  114. else:
  115. return expr
  116. if len(expr.free_symbols) < 2:
  117. return expr
  118. # don't destroy a Mul since much of the work may already be done
  119. if expr.is_Mul:
  120. args = list(expr.args)
  121. changed = False
  122. for i, a in enumerate(args):
  123. args[i] = separatevars(a, force)
  124. changed = changed or args[i] != a
  125. if changed:
  126. expr = expr.func(*args)
  127. return expr
  128. # get a Pow ready for expansion
  129. if expr.is_Pow and expr.base != S.Exp1:
  130. expr = Pow(separatevars(expr.base, force=force), expr.exp)
  131. # First try other expansion methods
  132. expr = expr.expand(mul=False, multinomial=False, force=force)
  133. _expr, reps = posify(expr) if force else (expr, {})
  134. expr = factor(_expr).subs(reps)
  135. if not expr.is_Add:
  136. return expr
  137. # Find any common coefficients to pull out
  138. args = list(expr.args)
  139. commonc = args[0].args_cnc(cset=True, warn=False)[0]
  140. for i in args[1:]:
  141. commonc &= i.args_cnc(cset=True, warn=False)[0]
  142. commonc = Mul(*commonc)
  143. commonc = commonc.as_coeff_Mul()[1] # ignore constants
  144. commonc_set = commonc.args_cnc(cset=True, warn=False)[0]
  145. # remove them
  146. for i, a in enumerate(args):
  147. c, nc = a.args_cnc(cset=True, warn=False)
  148. c = c - commonc_set
  149. args[i] = Mul(*c)*Mul(*nc)
  150. nonsepar = Add(*args)
  151. if len(nonsepar.free_symbols) > 1:
  152. _expr = nonsepar
  153. _expr, reps = posify(_expr) if force else (_expr, {})
  154. _expr = (factor(_expr)).subs(reps)
  155. if not _expr.is_Add:
  156. nonsepar = _expr
  157. return commonc*nonsepar
  158. def _separatevars_dict(expr, symbols):
  159. if symbols:
  160. if not all(t.is_Atom for t in symbols):
  161. raise ValueError("symbols must be Atoms.")
  162. symbols = list(symbols)
  163. elif symbols is None:
  164. return {'coeff': expr}
  165. else:
  166. symbols = list(expr.free_symbols)
  167. if not symbols:
  168. return None
  169. ret = {i: [] for i in symbols + ['coeff']}
  170. for i in Mul.make_args(expr):
  171. expsym = i.free_symbols
  172. intersection = set(symbols).intersection(expsym)
  173. if len(intersection) > 1:
  174. return None
  175. if len(intersection) == 0:
  176. # There are no symbols, so it is part of the coefficient
  177. ret['coeff'].append(i)
  178. else:
  179. ret[intersection.pop()].append(i)
  180. # rebuild
  181. for k, v in ret.items():
  182. ret[k] = Mul(*v)
  183. return ret
  184. def posify(eq):
  185. """Return ``eq`` (with generic symbols made positive) and a
  186. dictionary containing the mapping between the old and new
  187. symbols.
  188. Explanation
  189. ===========
  190. Any symbol that has positive=None will be replaced with a positive dummy
  191. symbol having the same name. This replacement will allow more symbolic
  192. processing of expressions, especially those involving powers and
  193. logarithms.
  194. A dictionary that can be sent to subs to restore ``eq`` to its original
  195. symbols is also returned.
  196. >>> from sympy import posify, Symbol, log, solve
  197. >>> from sympy.abc import x
  198. >>> posify(x + Symbol('p', positive=True) + Symbol('n', negative=True))
  199. (_x + n + p, {_x: x})
  200. >>> eq = 1/x
  201. >>> log(eq).expand()
  202. log(1/x)
  203. >>> log(posify(eq)[0]).expand()
  204. -log(_x)
  205. >>> p, rep = posify(eq)
  206. >>> log(p).expand().subs(rep)
  207. -log(x)
  208. It is possible to apply the same transformations to an iterable
  209. of expressions:
  210. >>> eq = x**2 - 4
  211. >>> solve(eq, x)
  212. [-2, 2]
  213. >>> eq_x, reps = posify([eq, x]); eq_x
  214. [_x**2 - 4, _x]
  215. >>> solve(*eq_x)
  216. [2]
  217. """
  218. eq = sympify(eq)
  219. if not isinstance(eq, Basic) and iterable(eq):
  220. f = type(eq)
  221. eq = list(eq)
  222. syms = set()
  223. for e in eq:
  224. syms = syms.union(e.atoms(Symbol))
  225. reps = {}
  226. for s in syms:
  227. reps.update({v: k for k, v in posify(s)[1].items()})
  228. for i, e in enumerate(eq):
  229. eq[i] = e.subs(reps)
  230. return f(eq), {r: s for s, r in reps.items()}
  231. reps = {s: Dummy(s.name, positive=True, **s.assumptions0)
  232. for s in eq.free_symbols if s.is_positive is None}
  233. eq = eq.subs(reps)
  234. return eq, {r: s for s, r in reps.items()}
  235. def hypersimp(f, k):
  236. """Given combinatorial term f(k) simplify its consecutive term ratio
  237. i.e. f(k+1)/f(k). The input term can be composed of functions and
  238. integer sequences which have equivalent representation in terms
  239. of gamma special function.
  240. Explanation
  241. ===========
  242. The algorithm performs three basic steps:
  243. 1. Rewrite all functions in terms of gamma, if possible.
  244. 2. Rewrite all occurrences of gamma in terms of products
  245. of gamma and rising factorial with integer, absolute
  246. constant exponent.
  247. 3. Perform simplification of nested fractions, powers
  248. and if the resulting expression is a quotient of
  249. polynomials, reduce their total degree.
  250. If f(k) is hypergeometric then as result we arrive with a
  251. quotient of polynomials of minimal degree. Otherwise None
  252. is returned.
  253. For more information on the implemented algorithm refer to:
  254. 1. W. Koepf, Algorithms for m-fold Hypergeometric Summation,
  255. Journal of Symbolic Computation (1995) 20, 399-417
  256. """
  257. f = sympify(f)
  258. g = f.subs(k, k + 1) / f
  259. g = g.rewrite(gamma)
  260. if g.has(Piecewise):
  261. g = piecewise_fold(g)
  262. g = g.args[-1][0]
  263. g = expand_func(g)
  264. g = powsimp(g, deep=True, combine='exp')
  265. if g.is_rational_function(k):
  266. return simplify(g, ratio=S.Infinity)
  267. else:
  268. return None
  269. def hypersimilar(f, g, k):
  270. """
  271. Returns True if ``f`` and ``g`` are hyper-similar.
  272. Explanation
  273. ===========
  274. Similarity in hypergeometric sense means that a quotient of
  275. f(k) and g(k) is a rational function in ``k``. This procedure
  276. is useful in solving recurrence relations.
  277. For more information see hypersimp().
  278. """
  279. f, g = list(map(sympify, (f, g)))
  280. h = (f/g).rewrite(gamma)
  281. h = h.expand(func=True, basic=False)
  282. return h.is_rational_function(k)
  283. def signsimp(expr, evaluate=None):
  284. """Make all Add sub-expressions canonical wrt sign.
  285. Explanation
  286. ===========
  287. If an Add subexpression, ``a``, can have a sign extracted,
  288. as determined by could_extract_minus_sign, it is replaced
  289. with Mul(-1, a, evaluate=False). This allows signs to be
  290. extracted from powers and products.
  291. Examples
  292. ========
  293. >>> from sympy import signsimp, exp, symbols
  294. >>> from sympy.abc import x, y
  295. >>> i = symbols('i', odd=True)
  296. >>> n = -1 + 1/x
  297. >>> n/x/(-n)**2 - 1/n/x
  298. (-1 + 1/x)/(x*(1 - 1/x)**2) - 1/(x*(-1 + 1/x))
  299. >>> signsimp(_)
  300. 0
  301. >>> x*n + x*-n
  302. x*(-1 + 1/x) + x*(1 - 1/x)
  303. >>> signsimp(_)
  304. 0
  305. Since powers automatically handle leading signs
  306. >>> (-2)**i
  307. -2**i
  308. signsimp can be used to put the base of a power with an integer
  309. exponent into canonical form:
  310. >>> n**i
  311. (-1 + 1/x)**i
  312. By default, signsimp does not leave behind any hollow simplification:
  313. if making an Add canonical wrt sign didn't change the expression, the
  314. original Add is restored. If this is not desired then the keyword
  315. ``evaluate`` can be set to False:
  316. >>> e = exp(y - x)
  317. >>> signsimp(e) == e
  318. True
  319. >>> signsimp(e, evaluate=False)
  320. exp(-(x - y))
  321. """
  322. if evaluate is None:
  323. evaluate = global_parameters.evaluate
  324. expr = sympify(expr)
  325. if not isinstance(expr, (Expr, Relational)) or expr.is_Atom:
  326. return expr
  327. # get rid of an pre-existing unevaluation regarding sign
  328. e = expr.replace(lambda x: x.is_Mul and -(-x) != x, lambda x: -(-x))
  329. e = sub_post(sub_pre(e))
  330. if not isinstance(e, (Expr, Relational)) or e.is_Atom:
  331. return e
  332. if e.is_Add:
  333. rv = e.func(*[signsimp(a) for a in e.args])
  334. if not evaluate and isinstance(rv, Add
  335. ) and rv.could_extract_minus_sign():
  336. return Mul(S.NegativeOne, -rv, evaluate=False)
  337. return rv
  338. if evaluate:
  339. e = e.replace(lambda x: x.is_Mul and -(-x) != x, lambda x: -(-x))
  340. return e
  341. @overload
  342. def simplify(expr: Expr, **kwargs) -> Expr: ...
  343. @overload
  344. def simplify(expr: Boolean, **kwargs) -> Boolean: ...
  345. @overload
  346. def simplify(expr: Set, **kwargs) -> Set: ...
  347. @overload
  348. def simplify(expr: Basic, **kwargs) -> Basic: ...
  349. def simplify(expr, ratio=1.7, measure=count_ops, rational=False, inverse=False, doit=True, **kwargs):
  350. """Simplifies the given expression.
  351. Explanation
  352. ===========
  353. Simplification is not a well defined term and the exact strategies
  354. this function tries can change in the future versions of SymPy. If
  355. your algorithm relies on "simplification" (whatever it is), try to
  356. determine what you need exactly - is it powsimp()?, radsimp()?,
  357. together()?, logcombine()?, or something else? And use this particular
  358. function directly, because those are well defined and thus your algorithm
  359. will be robust.
  360. Nonetheless, especially for interactive use, or when you do not know
  361. anything about the structure of the expression, simplify() tries to apply
  362. intelligent heuristics to make the input expression "simpler". For
  363. example:
  364. >>> from sympy import simplify, cos, sin
  365. >>> from sympy.abc import x, y
  366. >>> a = (x + x**2)/(x*sin(y)**2 + x*cos(y)**2)
  367. >>> a
  368. (x**2 + x)/(x*sin(y)**2 + x*cos(y)**2)
  369. >>> simplify(a)
  370. x + 1
  371. Note that we could have obtained the same result by using specific
  372. simplification functions:
  373. >>> from sympy import trigsimp, cancel
  374. >>> trigsimp(a)
  375. (x**2 + x)/x
  376. >>> cancel(_)
  377. x + 1
  378. In some cases, applying :func:`simplify` may actually result in some more
  379. complicated expression. The default ``ratio=1.7`` prevents more extreme
  380. cases: if (result length)/(input length) > ratio, then input is returned
  381. unmodified. The ``measure`` parameter lets you specify the function used
  382. to determine how complex an expression is. The function should take a
  383. single argument as an expression and return a number such that if
  384. expression ``a`` is more complex than expression ``b``, then
  385. ``measure(a) > measure(b)``. The default measure function is
  386. :func:`~.count_ops`, which returns the total number of operations in the
  387. expression.
  388. For example, if ``ratio=1``, ``simplify`` output cannot be longer
  389. than input.
  390. ::
  391. >>> from sympy import sqrt, simplify, count_ops, oo
  392. >>> root = 1/(sqrt(2)+3)
  393. Since ``simplify(root)`` would result in a slightly longer expression,
  394. root is returned unchanged instead::
  395. >>> simplify(root, ratio=1) == root
  396. True
  397. If ``ratio=oo``, simplify will be applied anyway::
  398. >>> count_ops(simplify(root, ratio=oo)) > count_ops(root)
  399. True
  400. Note that the shortest expression is not necessary the simplest, so
  401. setting ``ratio`` to 1 may not be a good idea.
  402. Heuristically, the default value ``ratio=1.7`` seems like a reasonable
  403. choice.
  404. You can easily define your own measure function based on what you feel
  405. should represent the "size" or "complexity" of the input expression. Note
  406. that some choices, such as ``lambda expr: len(str(expr))`` may appear to be
  407. good metrics, but have other problems (in this case, the measure function
  408. may slow down simplify too much for very large expressions). If you do not
  409. know what a good metric would be, the default, ``count_ops``, is a good
  410. one.
  411. For example:
  412. >>> from sympy import symbols, log
  413. >>> a, b = symbols('a b', positive=True)
  414. >>> g = log(a) + log(b) + log(a)*log(1/b)
  415. >>> h = simplify(g)
  416. >>> h
  417. log(a*b**(1 - log(a)))
  418. >>> count_ops(g)
  419. 8
  420. >>> count_ops(h)
  421. 5
  422. So you can see that ``h`` is simpler than ``g`` using the count_ops metric.
  423. However, we may not like how ``simplify`` (in this case, using
  424. ``logcombine``) has created the ``b**(log(1/a) + 1)`` term. A simple way
  425. to reduce this would be to give more weight to powers as operations in
  426. ``count_ops``. We can do this by using the ``visual=True`` option:
  427. >>> print(count_ops(g, visual=True))
  428. 2*ADD + DIV + 4*LOG + MUL
  429. >>> print(count_ops(h, visual=True))
  430. 2*LOG + MUL + POW + SUB
  431. >>> from sympy import Symbol, S
  432. >>> def my_measure(expr):
  433. ... POW = Symbol('POW')
  434. ... # Discourage powers by giving POW a weight of 10
  435. ... count = count_ops(expr, visual=True).subs(POW, 10)
  436. ... # Every other operation gets a weight of 1 (the default)
  437. ... count = count.replace(Symbol, type(S.One))
  438. ... return count
  439. >>> my_measure(g)
  440. 8
  441. >>> my_measure(h)
  442. 14
  443. >>> 15./8 > 1.7 # 1.7 is the default ratio
  444. True
  445. >>> simplify(g, measure=my_measure)
  446. -log(a)*log(b) + log(a) + log(b)
  447. Note that because ``simplify()`` internally tries many different
  448. simplification strategies and then compares them using the measure
  449. function, we get a completely different result that is still different
  450. from the input expression by doing this.
  451. If ``rational=True``, Floats will be recast as Rationals before simplification.
  452. If ``rational=None``, Floats will be recast as Rationals but the result will
  453. be recast as Floats. If rational=False(default) then nothing will be done
  454. to the Floats.
  455. If ``inverse=True``, it will be assumed that a composition of inverse
  456. functions, such as sin and asin, can be cancelled in any order.
  457. For example, ``asin(sin(x))`` will yield ``x`` without checking whether
  458. x belongs to the set where this relation is true. The default is
  459. False.
  460. Note that ``simplify()`` automatically calls ``doit()`` on the final
  461. expression. You can avoid this behavior by passing ``doit=False`` as
  462. an argument.
  463. Also, it should be noted that simplifying a boolean expression is not
  464. well defined. If the expression prefers automatic evaluation (such as
  465. :obj:`~.Eq()` or :obj:`~.Or()`), simplification will return ``True`` or
  466. ``False`` if truth value can be determined. If the expression is not
  467. evaluated by default (such as :obj:`~.Predicate()`), simplification will
  468. not reduce it and you should use :func:`~.refine` or :func:`~.ask`
  469. function. This inconsistency will be resolved in future version.
  470. See Also
  471. ========
  472. sympy.assumptions.refine.refine : Simplification using assumptions.
  473. sympy.assumptions.ask.ask : Query for boolean expressions using assumptions.
  474. """
  475. def shorter(*choices):
  476. """
  477. Return the choice that has the fewest ops. In case of a tie,
  478. the expression listed first is selected.
  479. """
  480. if not has_variety(choices):
  481. return choices[0]
  482. return min(choices, key=measure)
  483. def done(e):
  484. rv = e.doit() if doit else e
  485. return shorter(rv, collect_abs(rv))
  486. expr = sympify(expr, rational=rational)
  487. kwargs = {
  488. "ratio": kwargs.get('ratio', ratio),
  489. "measure": kwargs.get('measure', measure),
  490. "rational": kwargs.get('rational', rational),
  491. "inverse": kwargs.get('inverse', inverse),
  492. "doit": kwargs.get('doit', doit)}
  493. # no routine for Expr needs to check for is_zero
  494. if isinstance(expr, Expr) and expr.is_zero:
  495. return S.Zero if not expr.is_Number else expr
  496. _eval_simplify = getattr(expr, '_eval_simplify', None)
  497. if _eval_simplify is not None:
  498. return _eval_simplify(**kwargs)
  499. original_expr = expr = collect_abs(signsimp(expr))
  500. if not isinstance(expr, Basic) or not expr.args: # XXX: temporary hack
  501. return expr
  502. if inverse and expr.has(Function):
  503. expr = inversecombine(expr)
  504. if not expr.args: # simplified to atomic
  505. return expr
  506. # do deep simplification
  507. handled = Add, Mul, Pow, ExpBase
  508. expr = expr.replace(
  509. # here, checking for x.args is not enough because Basic has
  510. # args but Basic does not always play well with replace, e.g.
  511. # when simultaneous is True found expressions will be masked
  512. # off with a Dummy but not all Basic objects in an expression
  513. # can be replaced with a Dummy
  514. lambda x: isinstance(x, Expr) and x.args and not isinstance(
  515. x, handled),
  516. lambda x: x.func(*[simplify(i, **kwargs) for i in x.args]),
  517. simultaneous=False)
  518. if not isinstance(expr, handled):
  519. return done(expr)
  520. if not expr.is_commutative:
  521. expr = nc_simplify(expr)
  522. # TODO: Apply different strategies, considering expression pattern:
  523. # is it a purely rational function? Is there any trigonometric function?...
  524. # See also https://github.com/sympy/sympy/pull/185.
  525. # rationalize Floats
  526. floats = False
  527. if rational is not False and expr.has(Float):
  528. floats = True
  529. expr = nsimplify(expr, rational=True)
  530. expr = _bottom_up(expr, lambda w: getattr(w, 'normal', lambda: w)())
  531. expr = Mul(*powsimp(expr).as_content_primitive())
  532. _e = cancel(expr)
  533. expr1 = shorter(_e, _mexpand(_e).cancel()) # issue 6829
  534. expr2 = shorter(together(expr, deep=True), together(expr1, deep=True))
  535. if ratio is S.Infinity:
  536. expr = expr2
  537. else:
  538. expr = shorter(expr2, expr1, expr)
  539. if not isinstance(expr, Basic): # XXX: temporary hack
  540. return expr
  541. expr = factor_terms(expr, sign=False)
  542. # must come before `Piecewise` since this introduces more `Piecewise` terms
  543. if expr.has(sign):
  544. expr = expr.rewrite(Abs)
  545. # Deal with Piecewise separately to avoid recursive growth of expressions
  546. if expr.has(Piecewise):
  547. # Fold into a single Piecewise
  548. expr = piecewise_fold(expr)
  549. # Apply doit, if doit=True
  550. expr = done(expr)
  551. # Still a Piecewise?
  552. if expr.has(Piecewise):
  553. # Fold into a single Piecewise, in case doit lead to some
  554. # expressions being Piecewise
  555. expr = piecewise_fold(expr)
  556. # kroneckersimp also affects Piecewise
  557. if expr.has(KroneckerDelta):
  558. expr = kroneckersimp(expr)
  559. # Still a Piecewise?
  560. if expr.has(Piecewise):
  561. # Do not apply doit on the segments as it has already
  562. # been done above, but simplify
  563. expr = piecewise_simplify(expr, deep=True, doit=False)
  564. # Still a Piecewise?
  565. if expr.has(Piecewise):
  566. # Try factor common terms
  567. expr = shorter(expr, factor_terms(expr))
  568. # As all expressions have been simplified above with the
  569. # complete simplify, nothing more needs to be done here
  570. return expr
  571. # hyperexpand automatically only works on hypergeometric terms
  572. # Do this after the Piecewise part to avoid recursive expansion
  573. expr = hyperexpand(expr)
  574. if expr.has(KroneckerDelta):
  575. expr = kroneckersimp(expr)
  576. if expr.has(BesselBase):
  577. expr = besselsimp(expr)
  578. if expr.has(TrigonometricFunction, HyperbolicFunction):
  579. expr = trigsimp(expr, deep=True)
  580. if expr.has(log):
  581. expr = shorter(expand_log(expr, deep=True), logcombine(expr))
  582. if expr.has(CombinatorialFunction, gamma):
  583. # expression with gamma functions or non-integer arguments is
  584. # automatically passed to gammasimp
  585. expr = combsimp(expr)
  586. if expr.has(Sum):
  587. expr = sum_simplify(expr, **kwargs)
  588. if expr.has(Integral):
  589. expr = expr.xreplace({
  590. i: factor_terms(i) for i in expr.atoms(Integral)})
  591. if expr.has(Product):
  592. expr = product_simplify(expr, **kwargs)
  593. from sympy.physics.units import Quantity
  594. if expr.has(Quantity):
  595. from sympy.physics.units.util import quantity_simplify
  596. expr = quantity_simplify(expr)
  597. short = shorter(powsimp(expr, combine='exp', deep=True), powsimp(expr), expr)
  598. short = shorter(short, cancel(short))
  599. short = shorter(short, factor_terms(short), expand_power_exp(expand_mul(short)))
  600. if short.has(TrigonometricFunction, HyperbolicFunction, ExpBase, exp):
  601. short = exptrigsimp(short)
  602. # get rid of hollow 2-arg Mul factorization
  603. hollow_mul = Transform(
  604. lambda x: Mul(*x.args),
  605. lambda x:
  606. x.is_Mul and
  607. len(x.args) == 2 and
  608. x.args[0].is_Number and
  609. x.args[1].is_Add and
  610. x.is_commutative)
  611. expr = short.xreplace(hollow_mul)
  612. numer, denom = expr.as_numer_denom()
  613. if denom.is_Add:
  614. n, d = fraction(radsimp(1/denom, symbolic=False, max_terms=1))
  615. if n is not S.One:
  616. expr = (numer*n).expand()/d
  617. if expr.could_extract_minus_sign():
  618. n, d = fraction(expr)
  619. if d != 0:
  620. expr = signsimp(-n/(-d))
  621. if measure(expr) > ratio*measure(original_expr):
  622. expr = original_expr
  623. # restore floats
  624. if floats and rational is None:
  625. expr = nfloat(expr, exponent=False)
  626. return done(expr)
  627. def sum_simplify(s, **kwargs):
  628. """Main function for Sum simplification"""
  629. if not isinstance(s, Add):
  630. s = s.xreplace({a: sum_simplify(a, **kwargs)
  631. for a in s.atoms(Add) if a.has(Sum)})
  632. s = expand(s)
  633. if not isinstance(s, Add):
  634. return s
  635. terms = s.args
  636. s_t = [] # Sum Terms
  637. o_t = [] # Other Terms
  638. for term in terms:
  639. sum_terms, other = sift(Mul.make_args(term),
  640. lambda i: isinstance(i, Sum), binary=True)
  641. if not sum_terms:
  642. o_t.append(term)
  643. continue
  644. other = [Mul(*other)]
  645. s_t.append(Mul(*(other + [s._eval_simplify(**kwargs) for s in sum_terms])))
  646. result = Add(sum_combine(s_t), *o_t)
  647. return result
  648. def sum_combine(s_t):
  649. """Helper function for Sum simplification
  650. Attempts to simplify a list of sums, by combining limits / sum function's
  651. returns the simplified sum
  652. """
  653. used = [False] * len(s_t)
  654. for method in range(2):
  655. for i, s_term1 in enumerate(s_t):
  656. if not used[i]:
  657. for j, s_term2 in enumerate(s_t):
  658. if not used[j] and i != j:
  659. temp = sum_add(s_term1, s_term2, method)
  660. if isinstance(temp, (Sum, Mul)):
  661. s_t[i] = temp
  662. s_term1 = s_t[i]
  663. used[j] = True
  664. result = S.Zero
  665. for i, s_term in enumerate(s_t):
  666. if not used[i]:
  667. result = Add(result, s_term)
  668. return result
  669. def factor_sum(self, limits=None, radical=False, clear=False, fraction=False, sign=True):
  670. """Return Sum with constant factors extracted.
  671. If ``limits`` is specified then ``self`` is the summand; the other
  672. keywords are passed to ``factor_terms``.
  673. Examples
  674. ========
  675. >>> from sympy import Sum
  676. >>> from sympy.abc import x, y
  677. >>> from sympy.simplify.simplify import factor_sum
  678. >>> s = Sum(x*y, (x, 1, 3))
  679. >>> factor_sum(s)
  680. y*Sum(x, (x, 1, 3))
  681. >>> factor_sum(s.function, s.limits)
  682. y*Sum(x, (x, 1, 3))
  683. """
  684. # XXX deprecate in favor of direct call to factor_terms
  685. kwargs = {"radical": radical, "clear": clear,
  686. "fraction": fraction, "sign": sign}
  687. expr = Sum(self, *limits) if limits else self
  688. return factor_terms(expr, **kwargs)
  689. def sum_add(self, other, method=0):
  690. """Helper function for Sum simplification"""
  691. #we know this is something in terms of a constant * a sum
  692. #so we temporarily put the constants inside for simplification
  693. #then simplify the result
  694. def __refactor(val):
  695. args = Mul.make_args(val)
  696. sumv = next(x for x in args if isinstance(x, Sum))
  697. constant = Mul(*[x for x in args if x != sumv])
  698. return Sum(constant * sumv.function, *sumv.limits)
  699. if isinstance(self, Mul):
  700. rself = __refactor(self)
  701. else:
  702. rself = self
  703. if isinstance(other, Mul):
  704. rother = __refactor(other)
  705. else:
  706. rother = other
  707. if type(rself) is type(rother):
  708. if method == 0:
  709. if rself.limits == rother.limits:
  710. return factor_sum(Sum(rself.function + rother.function, *rself.limits))
  711. elif method == 1:
  712. if simplify(rself.function - rother.function) == 0:
  713. if len(rself.limits) == len(rother.limits) == 1:
  714. i = rself.limits[0][0]
  715. x1 = rself.limits[0][1]
  716. y1 = rself.limits[0][2]
  717. j = rother.limits[0][0]
  718. x2 = rother.limits[0][1]
  719. y2 = rother.limits[0][2]
  720. if i == j:
  721. if x2 == y1 + 1:
  722. return factor_sum(Sum(rself.function, (i, x1, y2)))
  723. elif x1 == y2 + 1:
  724. return factor_sum(Sum(rself.function, (i, x2, y1)))
  725. return Add(self, other)
  726. def product_simplify(s, **kwargs):
  727. """Main function for Product simplification"""
  728. terms = Mul.make_args(s)
  729. p_t = [] # Product Terms
  730. o_t = [] # Other Terms
  731. deep = kwargs.get('deep', True)
  732. for term in terms:
  733. if isinstance(term, Product):
  734. if deep:
  735. p_t.append(Product(term.function.simplify(**kwargs),
  736. *term.limits))
  737. else:
  738. p_t.append(term)
  739. else:
  740. o_t.append(term)
  741. used = [False] * len(p_t)
  742. for method in range(2):
  743. for i, p_term1 in enumerate(p_t):
  744. if not used[i]:
  745. for j, p_term2 in enumerate(p_t):
  746. if not used[j] and i != j:
  747. tmp_prod = product_mul(p_term1, p_term2, method)
  748. if isinstance(tmp_prod, Product):
  749. p_t[i] = tmp_prod
  750. used[j] = True
  751. result = Mul(*o_t)
  752. for i, p_term in enumerate(p_t):
  753. if not used[i]:
  754. result = Mul(result, p_term)
  755. return result
  756. def product_mul(self, other, method=0):
  757. """Helper function for Product simplification"""
  758. if type(self) is type(other):
  759. if method == 0:
  760. if self.limits == other.limits:
  761. return Product(self.function * other.function, *self.limits)
  762. elif method == 1:
  763. if simplify(self.function - other.function) == 0:
  764. if len(self.limits) == len(other.limits) == 1:
  765. i = self.limits[0][0]
  766. x1 = self.limits[0][1]
  767. y1 = self.limits[0][2]
  768. j = other.limits[0][0]
  769. x2 = other.limits[0][1]
  770. y2 = other.limits[0][2]
  771. if i == j:
  772. if x2 == y1 + 1:
  773. return Product(self.function, (i, x1, y2))
  774. elif x1 == y2 + 1:
  775. return Product(self.function, (i, x2, y1))
  776. return Mul(self, other)
  777. def _nthroot_solve(p, n, prec):
  778. """
  779. helper function for ``nthroot``
  780. It denests ``p**Rational(1, n)`` using its minimal polynomial
  781. """
  782. from sympy.solvers import solve
  783. while n % 2 == 0:
  784. p = sqrtdenest(sqrt(p))
  785. n = n // 2
  786. if n == 1:
  787. return p
  788. pn = p**Rational(1, n)
  789. x = Symbol('x')
  790. f = _minimal_polynomial_sq(p, n, x)
  791. if f is None:
  792. return None
  793. sols = solve(f, x)
  794. for sol in sols:
  795. if abs(sol - pn).n() < 1./10**prec:
  796. sol = sqrtdenest(sol)
  797. if _mexpand(sol**n) == p:
  798. return sol
  799. def logcombine(expr, force=False):
  800. """
  801. Takes logarithms and combines them using the following rules:
  802. - log(x) + log(y) == log(x*y) if both are positive
  803. - a*log(x) == log(x**a) if x is positive and a is real
  804. If ``force`` is ``True`` then the assumptions above will be assumed to hold if
  805. there is no assumption already in place on a quantity. For example, if
  806. ``a`` is imaginary or the argument negative, force will not perform a
  807. combination but if ``a`` is a symbol with no assumptions the change will
  808. take place.
  809. Examples
  810. ========
  811. >>> from sympy import Symbol, symbols, log, logcombine, I
  812. >>> from sympy.abc import a, x, y, z
  813. >>> logcombine(a*log(x) + log(y) - log(z))
  814. a*log(x) + log(y) - log(z)
  815. >>> logcombine(a*log(x) + log(y) - log(z), force=True)
  816. log(x**a*y/z)
  817. >>> x,y,z = symbols('x,y,z', positive=True)
  818. >>> a = Symbol('a', real=True)
  819. >>> logcombine(a*log(x) + log(y) - log(z))
  820. log(x**a*y/z)
  821. The transformation is limited to factors and/or terms that
  822. contain logs, so the result depends on the initial state of
  823. expansion:
  824. >>> eq = (2 + 3*I)*log(x)
  825. >>> logcombine(eq, force=True) == eq
  826. True
  827. >>> logcombine(eq.expand(), force=True)
  828. log(x**2) + I*log(x**3)
  829. See Also
  830. ========
  831. posify: replace all symbols with symbols having positive assumptions
  832. sympy.core.function.expand_log: expand the logarithms of products
  833. and powers; the opposite of logcombine
  834. """
  835. def f(rv):
  836. if not (rv.is_Add or rv.is_Mul):
  837. return rv
  838. def gooda(a):
  839. # bool to tell whether the leading ``a`` in ``a*log(x)``
  840. # could appear as log(x**a)
  841. return (a is not S.NegativeOne and # -1 *could* go, but we disallow
  842. (a.is_extended_real or force and a.is_extended_real is not False))
  843. def goodlog(l):
  844. # bool to tell whether log ``l``'s argument can combine with others
  845. a = l.args[0]
  846. return a.is_positive or force and a.is_nonpositive is not False
  847. other = []
  848. logs = []
  849. log1 = defaultdict(list)
  850. for a in Add.make_args(rv):
  851. if isinstance(a, log) and goodlog(a):
  852. log1[()].append(([], a))
  853. elif not a.is_Mul:
  854. other.append(a)
  855. else:
  856. ot = []
  857. co = []
  858. lo = []
  859. for ai in a.args:
  860. if ai.is_Rational and ai < 0:
  861. ot.append(S.NegativeOne)
  862. co.append(-ai)
  863. elif isinstance(ai, log) and goodlog(ai):
  864. lo.append(ai)
  865. elif gooda(ai):
  866. co.append(ai)
  867. else:
  868. ot.append(ai)
  869. if len(lo) > 1:
  870. logs.append((ot, co, lo))
  871. elif lo:
  872. log1[tuple(ot)].append((co, lo[0]))
  873. else:
  874. other.append(a)
  875. # if there is only one log in other, put it with the
  876. # good logs
  877. if len(other) == 1 and isinstance(other[0], log):
  878. log1[()].append(([], other.pop()))
  879. # if there is only one log at each coefficient and none have
  880. # an exponent to place inside the log then there is nothing to do
  881. if not logs and all(len(log1[k]) == 1 and log1[k][0] == [] for k in log1):
  882. return rv
  883. # collapse multi-logs as far as possible in a canonical way
  884. # TODO: see if x*log(a)+x*log(a)*log(b) -> x*log(a)*(1+log(b))?
  885. # -- in this case, it's unambiguous, but if it were were a log(c) in
  886. # each term then it's arbitrary whether they are grouped by log(a) or
  887. # by log(c). So for now, just leave this alone; it's probably better to
  888. # let the user decide
  889. for o, e, l in logs:
  890. l = list(ordered(l))
  891. e = log(l.pop(0).args[0]**Mul(*e))
  892. while l:
  893. li = l.pop(0)
  894. e = log(li.args[0]**e)
  895. c, l = Mul(*o), e
  896. if isinstance(l, log): # it should be, but check to be sure
  897. log1[(c,)].append(([], l))
  898. else:
  899. other.append(c*l)
  900. # logs that have the same coefficient can multiply
  901. for k in list(log1.keys()):
  902. log1[Mul(*k)] = log(logcombine(Mul(*[
  903. l.args[0]**Mul(*c) for c, l in log1.pop(k)]),
  904. force=force), evaluate=False)
  905. # logs that have oppositely signed coefficients can divide
  906. for k in ordered(list(log1.keys())):
  907. if k not in log1: # already popped as -k
  908. continue
  909. if -k in log1:
  910. # figure out which has the minus sign; the one with
  911. # more op counts should be the one
  912. num, den = k, -k
  913. if num.count_ops() > den.count_ops():
  914. num, den = den, num
  915. other.append(
  916. num*log(log1.pop(num).args[0]/log1.pop(den).args[0],
  917. evaluate=False))
  918. else:
  919. other.append(k*log1.pop(k))
  920. return Add(*other)
  921. return _bottom_up(expr, f)
  922. def inversecombine(expr):
  923. """Simplify the composition of a function and its inverse.
  924. Explanation
  925. ===========
  926. No attention is paid to whether the inverse is a left inverse or a
  927. right inverse; thus, the result will in general not be equivalent
  928. to the original expression.
  929. Examples
  930. ========
  931. >>> from sympy.simplify.simplify import inversecombine
  932. >>> from sympy import asin, sin, log, exp
  933. >>> from sympy.abc import x
  934. >>> inversecombine(asin(sin(x)))
  935. x
  936. >>> inversecombine(2*log(exp(3*x)))
  937. 6*x
  938. """
  939. def f(rv):
  940. if isinstance(rv, log):
  941. if isinstance(rv.args[0], exp) or (rv.args[0].is_Pow and rv.args[0].base == S.Exp1):
  942. rv = rv.args[0].exp
  943. elif rv.is_Function and hasattr(rv, "inverse"):
  944. if (len(rv.args) == 1 and len(rv.args[0].args) == 1 and
  945. isinstance(rv.args[0], rv.inverse(argindex=1))):
  946. rv = rv.args[0].args[0]
  947. if rv.is_Pow and rv.base == S.Exp1:
  948. if isinstance(rv.exp, log):
  949. rv = rv.exp.args[0]
  950. return rv
  951. return _bottom_up(expr, f)
  952. def kroneckersimp(expr):
  953. """
  954. Simplify expressions with KroneckerDelta.
  955. The only simplification currently attempted is to identify multiplicative cancellation:
  956. Examples
  957. ========
  958. >>> from sympy import KroneckerDelta, kroneckersimp
  959. >>> from sympy.abc import i
  960. >>> kroneckersimp(1 + KroneckerDelta(0, i) * KroneckerDelta(1, i))
  961. 1
  962. """
  963. def args_cancel(args1, args2):
  964. for i1 in range(2):
  965. for i2 in range(2):
  966. a1 = args1[i1]
  967. a2 = args2[i2]
  968. a3 = args1[(i1 + 1) % 2]
  969. a4 = args2[(i2 + 1) % 2]
  970. if Eq(a1, a2) is S.true and Eq(a3, a4) is S.false:
  971. return True
  972. return False
  973. def cancel_kronecker_mul(m):
  974. args = m.args
  975. deltas = [a for a in args if isinstance(a, KroneckerDelta)]
  976. for delta1, delta2 in subsets(deltas, 2):
  977. args1 = delta1.args
  978. args2 = delta2.args
  979. if args_cancel(args1, args2):
  980. return S.Zero * m # In case of oo etc
  981. return m
  982. if not expr.has(KroneckerDelta):
  983. return expr
  984. if expr.has(Piecewise):
  985. expr = expr.rewrite(KroneckerDelta)
  986. newexpr = expr
  987. expr = None
  988. while newexpr != expr:
  989. expr = newexpr
  990. newexpr = expr.replace(lambda e: isinstance(e, Mul), cancel_kronecker_mul)
  991. return expr
  992. def besselsimp(expr):
  993. """
  994. Simplify bessel-type functions.
  995. Explanation
  996. ===========
  997. This routine tries to simplify bessel-type functions. Currently it only
  998. works on the Bessel J and I functions, however. It works by looking at all
  999. such functions in turn, and eliminating factors of "I" and "-1" (actually
  1000. their polar equivalents) in front of the argument. Then, functions of
  1001. half-integer order are rewritten using trigonometric functions and
  1002. functions of integer order (> 1) are rewritten using functions
  1003. of low order. Finally, if the expression was changed, compute
  1004. factorization of the result with factor().
  1005. >>> from sympy import besselj, besseli, besselsimp, polar_lift, I, S
  1006. >>> from sympy.abc import z, nu
  1007. >>> besselsimp(besselj(nu, z*polar_lift(-1)))
  1008. exp(I*pi*nu)*besselj(nu, z)
  1009. >>> besselsimp(besseli(nu, z*polar_lift(-I)))
  1010. exp(-I*pi*nu/2)*besselj(nu, z)
  1011. >>> besselsimp(besseli(S(-1)/2, z))
  1012. sqrt(2)*cosh(z)/(sqrt(pi)*sqrt(z))
  1013. >>> besselsimp(z*besseli(0, z) + z*(besseli(2, z))/2 + besseli(1, z))
  1014. 3*z*besseli(0, z)/2
  1015. """
  1016. # TODO
  1017. # - better algorithm?
  1018. # - simplify (cos(pi*b)*besselj(b,z) - besselj(-b,z))/sin(pi*b) ...
  1019. # - use contiguity relations?
  1020. def replacer(fro, to, factors):
  1021. factors = set(factors)
  1022. def repl(nu, z):
  1023. if factors.intersection(Mul.make_args(z)):
  1024. return to(nu, z)
  1025. return fro(nu, z)
  1026. return repl
  1027. def torewrite(fro, to):
  1028. def tofunc(nu, z):
  1029. return fro(nu, z).rewrite(to)
  1030. return tofunc
  1031. def tominus(fro):
  1032. def tofunc(nu, z):
  1033. return exp(I*pi*nu)*fro(nu, exp_polar(-I*pi)*z)
  1034. return tofunc
  1035. orig_expr = expr
  1036. ifactors = [I, exp_polar(I*pi/2), exp_polar(-I*pi/2)]
  1037. expr = expr.replace(
  1038. besselj, replacer(besselj,
  1039. torewrite(besselj, besseli), ifactors))
  1040. expr = expr.replace(
  1041. besseli, replacer(besseli,
  1042. torewrite(besseli, besselj), ifactors))
  1043. minusfactors = [-1, exp_polar(I*pi)]
  1044. expr = expr.replace(
  1045. besselj, replacer(besselj, tominus(besselj), minusfactors))
  1046. expr = expr.replace(
  1047. besseli, replacer(besseli, tominus(besseli), minusfactors))
  1048. z0 = Dummy('z')
  1049. def expander(fro):
  1050. def repl(nu, z):
  1051. if (nu % 1) == S.Half:
  1052. return simplify(trigsimp(unpolarify(
  1053. fro(nu, z0).rewrite(besselj).rewrite(jn).expand(
  1054. func=True)).subs(z0, z)))
  1055. elif nu.is_Integer and nu > 1:
  1056. return fro(nu, z).expand(func=True)
  1057. return fro(nu, z)
  1058. return repl
  1059. expr = expr.replace(besselj, expander(besselj))
  1060. expr = expr.replace(bessely, expander(bessely))
  1061. expr = expr.replace(besseli, expander(besseli))
  1062. expr = expr.replace(besselk, expander(besselk))
  1063. def _bessel_simp_recursion(expr):
  1064. def _use_recursion(bessel, expr):
  1065. while True:
  1066. bessels = expr.find(lambda x: isinstance(x, bessel))
  1067. try:
  1068. for ba in sorted(bessels, key=lambda x: re(x.args[0])):
  1069. a, x = ba.args
  1070. bap1 = bessel(a+1, x)
  1071. bap2 = bessel(a+2, x)
  1072. if expr.has(bap1) and expr.has(bap2):
  1073. expr = expr.subs(ba, 2*(a+1)/x*bap1 - bap2)
  1074. break
  1075. else:
  1076. return expr
  1077. except (ValueError, TypeError):
  1078. return expr
  1079. if expr.has(besselj):
  1080. expr = _use_recursion(besselj, expr)
  1081. if expr.has(bessely):
  1082. expr = _use_recursion(bessely, expr)
  1083. return expr
  1084. expr = _bessel_simp_recursion(expr)
  1085. if expr != orig_expr:
  1086. expr = expr.factor()
  1087. return expr
  1088. def nthroot(expr, n, max_len=4, prec=15):
  1089. """
  1090. Compute a real nth-root of a sum of surds.
  1091. Parameters
  1092. ==========
  1093. expr : sum of surds
  1094. n : integer
  1095. max_len : maximum number of surds passed as constants to ``nsimplify``
  1096. Algorithm
  1097. =========
  1098. First ``nsimplify`` is used to get a candidate root; if it is not a
  1099. root the minimal polynomial is computed; the answer is one of its
  1100. roots.
  1101. Examples
  1102. ========
  1103. >>> from sympy.simplify.simplify import nthroot
  1104. >>> from sympy import sqrt
  1105. >>> nthroot(90 + 34*sqrt(7), 3)
  1106. sqrt(7) + 3
  1107. """
  1108. expr = sympify(expr)
  1109. n = sympify(n)
  1110. p = expr**Rational(1, n)
  1111. if not n.is_integer:
  1112. return p
  1113. if not _is_sum_surds(expr):
  1114. return p
  1115. surds = []
  1116. coeff_muls = [x.as_coeff_Mul() for x in expr.args]
  1117. for x, y in coeff_muls:
  1118. if not x.is_rational:
  1119. return p
  1120. if y is S.One:
  1121. continue
  1122. if not (y.is_Pow and y.exp == S.Half and y.base.is_integer):
  1123. return p
  1124. surds.append(y)
  1125. surds.sort()
  1126. surds = surds[:max_len]
  1127. if expr < 0 and n % 2 == 1:
  1128. p = (-expr)**Rational(1, n)
  1129. a = nsimplify(p, constants=surds)
  1130. res = a if _mexpand(a**n) == _mexpand(-expr) else p
  1131. return -res
  1132. a = nsimplify(p, constants=surds)
  1133. if _mexpand(a) is not _mexpand(p) and _mexpand(a**n) == _mexpand(expr):
  1134. return _mexpand(a)
  1135. expr = _nthroot_solve(expr, n, prec)
  1136. if expr is None:
  1137. return p
  1138. return expr
  1139. def nsimplify(expr, constants=(), tolerance=None, full=False, rational=None,
  1140. rational_conversion='base10'):
  1141. """
  1142. Find a simple representation for a number or, if there are free symbols or
  1143. if ``rational=True``, then replace Floats with their Rational equivalents. If
  1144. no change is made and rational is not False then Floats will at least be
  1145. converted to Rationals.
  1146. Explanation
  1147. ===========
  1148. For numerical expressions, a simple formula that numerically matches the
  1149. given numerical expression is sought (and the input should be possible
  1150. to evalf to a precision of at least 30 digits).
  1151. Optionally, a list of (rationally independent) constants to
  1152. include in the formula may be given.
  1153. A lower tolerance may be set to find less exact matches. If no tolerance
  1154. is given then the least precise value will set the tolerance (e.g. Floats
  1155. default to 15 digits of precision, so would be tolerance=10**-15).
  1156. With ``full=True``, a more extensive search is performed
  1157. (this is useful to find simpler numbers when the tolerance
  1158. is set low).
  1159. When converting to rational, if rational_conversion='base10' (the default), then
  1160. convert floats to rationals using their base-10 (string) representation.
  1161. When rational_conversion='exact' it uses the exact, base-2 representation.
  1162. Examples
  1163. ========
  1164. >>> from sympy import nsimplify, sqrt, GoldenRatio, exp, I, pi
  1165. >>> nsimplify(4/(1+sqrt(5)), [GoldenRatio])
  1166. -2 + 2*GoldenRatio
  1167. >>> nsimplify((1/(exp(3*pi*I/5)+1)))
  1168. 1/2 - I*sqrt(sqrt(5)/10 + 1/4)
  1169. >>> nsimplify(I**I, [pi])
  1170. exp(-pi/2)
  1171. >>> nsimplify(pi, tolerance=0.01)
  1172. 22/7
  1173. >>> nsimplify(0.333333333333333, rational=True, rational_conversion='exact')
  1174. 6004799503160655/18014398509481984
  1175. >>> nsimplify(0.333333333333333, rational=True)
  1176. 1/3
  1177. See Also
  1178. ========
  1179. sympy.core.function.nfloat
  1180. """
  1181. try:
  1182. return sympify(as_int(expr))
  1183. except (TypeError, ValueError):
  1184. pass
  1185. expr = sympify(expr).xreplace({
  1186. Float('inf'): S.Infinity,
  1187. Float('-inf'): S.NegativeInfinity,
  1188. })
  1189. if expr is S.Infinity or expr is S.NegativeInfinity:
  1190. return expr
  1191. if rational or expr.free_symbols:
  1192. return _real_to_rational(expr, tolerance, rational_conversion)
  1193. # SymPy's default tolerance for Rationals is 15; other numbers may have
  1194. # lower tolerances set, so use them to pick the largest tolerance if None
  1195. # was given
  1196. if tolerance is None:
  1197. tolerance = 10**-min([15] +
  1198. [mpmath.libmp.libmpf.prec_to_dps(n._prec)
  1199. for n in expr.atoms(Float)])
  1200. # XXX should prec be set independent of tolerance or should it be computed
  1201. # from tolerance?
  1202. prec = 30
  1203. bprec = int(prec*3.33)
  1204. constants_dict = {}
  1205. for constant in constants:
  1206. constant = sympify(constant)
  1207. v = constant.evalf(prec)
  1208. if not v.is_Float:
  1209. raise ValueError("constants must be real-valued")
  1210. constants_dict[str(constant)] = v._to_mpmath(bprec)
  1211. exprval = expr.evalf(prec, chop=True)
  1212. re, im = exprval.as_real_imag()
  1213. # safety check to make sure that this evaluated to a number
  1214. if not (re.is_Number and im.is_Number):
  1215. return expr
  1216. def nsimplify_real(x):
  1217. orig = mpmath.mp.dps
  1218. xv = x._to_mpmath(bprec)
  1219. try:
  1220. # We'll be happy with low precision if a simple fraction
  1221. if not (tolerance or full):
  1222. mpmath.mp.dps = 15
  1223. rat = mpmath.pslq([xv, 1])
  1224. if rat is not None:
  1225. return Rational(-int(rat[1]), int(rat[0]))
  1226. mpmath.mp.dps = prec
  1227. newexpr = mpmath.identify(xv, constants=constants_dict,
  1228. tol=tolerance, full=full)
  1229. if not newexpr:
  1230. raise ValueError
  1231. if full:
  1232. newexpr = newexpr[0]
  1233. expr = sympify(newexpr)
  1234. if x and not expr: # don't let x become 0
  1235. raise ValueError
  1236. if expr.is_finite is False and xv not in [mpmath.inf, mpmath.ninf]:
  1237. raise ValueError
  1238. return expr
  1239. finally:
  1240. # even though there are returns above, this is executed
  1241. # before leaving
  1242. mpmath.mp.dps = orig
  1243. try:
  1244. if re:
  1245. re = nsimplify_real(re)
  1246. if im:
  1247. im = nsimplify_real(im)
  1248. except ValueError:
  1249. if rational is None:
  1250. return _real_to_rational(expr, rational_conversion=rational_conversion)
  1251. return expr
  1252. rv = re + im*S.ImaginaryUnit
  1253. # if there was a change or rational is explicitly not wanted
  1254. # return the value, else return the Rational representation
  1255. if rv != expr or rational is False:
  1256. return rv
  1257. return _real_to_rational(expr, rational_conversion=rational_conversion)
  1258. def _real_to_rational(expr, tolerance=None, rational_conversion='base10'):
  1259. """
  1260. Replace all reals in expr with rationals.
  1261. Examples
  1262. ========
  1263. >>> from sympy.simplify.simplify import _real_to_rational
  1264. >>> from sympy.abc import x
  1265. >>> _real_to_rational(.76 + .1*x**.5)
  1266. sqrt(x)/10 + 19/25
  1267. If rational_conversion='base10', this uses the base-10 string. If
  1268. rational_conversion='exact', the exact, base-2 representation is used.
  1269. >>> _real_to_rational(0.333333333333333, rational_conversion='exact')
  1270. 6004799503160655/18014398509481984
  1271. >>> _real_to_rational(0.333333333333333)
  1272. 1/3
  1273. """
  1274. expr = _sympify(expr)
  1275. inf = Float('inf')
  1276. p = expr
  1277. reps = {}
  1278. reduce_num = None
  1279. if tolerance is not None and tolerance < 1:
  1280. reduce_num = ceiling(1/tolerance)
  1281. for fl in p.atoms(Float):
  1282. key = fl
  1283. if reduce_num is not None:
  1284. r = Rational(fl).limit_denominator(reduce_num)
  1285. elif (tolerance is not None and tolerance >= 1 and
  1286. fl.is_Integer is False):
  1287. r = Rational(tolerance*round(fl/tolerance)
  1288. ).limit_denominator(int(tolerance))
  1289. else:
  1290. if rational_conversion == 'exact':
  1291. r = Rational(fl)
  1292. reps[key] = r
  1293. continue
  1294. elif rational_conversion != 'base10':
  1295. raise ValueError("rational_conversion must be 'base10' or 'exact'")
  1296. r = nsimplify(fl, rational=False)
  1297. # e.g. log(3).n() -> log(3) instead of a Rational
  1298. if fl and not r:
  1299. r = Rational(fl)
  1300. elif not r.is_Rational:
  1301. if fl in (inf, -inf):
  1302. r = S.ComplexInfinity
  1303. elif fl < 0:
  1304. fl = -fl
  1305. d = Pow(10, int(mpmath.log(fl)/mpmath.log(10)))
  1306. r = -Rational(str(fl/d))*d
  1307. elif fl > 0:
  1308. d = Pow(10, int(mpmath.log(fl)/mpmath.log(10)))
  1309. r = Rational(str(fl/d))*d
  1310. else:
  1311. r = S.Zero
  1312. reps[key] = r
  1313. return p.subs(reps, simultaneous=True)
  1314. def clear_coefficients(expr, rhs=S.Zero):
  1315. """Return `p, r` where `p` is the expression obtained when Rational
  1316. additive and multiplicative coefficients of `expr` have been stripped
  1317. away in a naive fashion (i.e. without simplification). The operations
  1318. needed to remove the coefficients will be applied to `rhs` and returned
  1319. as `r`.
  1320. Examples
  1321. ========
  1322. >>> from sympy.simplify.simplify import clear_coefficients
  1323. >>> from sympy.abc import x, y
  1324. >>> from sympy import Dummy
  1325. >>> expr = 4*y*(6*x + 3)
  1326. >>> clear_coefficients(expr - 2)
  1327. (y*(2*x + 1), 1/6)
  1328. When solving 2 or more expressions like `expr = a`,
  1329. `expr = b`, etc..., it is advantageous to provide a Dummy symbol
  1330. for `rhs` and simply replace it with `a`, `b`, etc... in `r`.
  1331. >>> rhs = Dummy('rhs')
  1332. >>> clear_coefficients(expr, rhs)
  1333. (y*(2*x + 1), _rhs/12)
  1334. >>> _[1].subs(rhs, 2)
  1335. 1/6
  1336. """
  1337. was = None
  1338. free = expr.free_symbols
  1339. if expr.is_Rational:
  1340. return (S.Zero, rhs - expr)
  1341. while expr and was != expr:
  1342. was = expr
  1343. m, expr = (
  1344. expr.as_content_primitive()
  1345. if free else
  1346. factor_terms(expr).as_coeff_Mul(rational=True))
  1347. rhs /= m
  1348. c, expr = expr.as_coeff_Add(rational=True)
  1349. rhs -= c
  1350. expr = signsimp(expr, evaluate = False)
  1351. if expr.could_extract_minus_sign():
  1352. expr = -expr
  1353. rhs = -rhs
  1354. return expr, rhs
  1355. def nc_simplify(expr, deep=True):
  1356. '''
  1357. Simplify a non-commutative expression composed of multiplication
  1358. and raising to a power by grouping repeated subterms into one power.
  1359. Priority is given to simplifications that give the fewest number
  1360. of arguments in the end (for example, in a*b*a*b*c*a*b*c simplifying
  1361. to (a*b)**2*c*a*b*c gives 5 arguments while a*b*(a*b*c)**2 has 3).
  1362. If ``expr`` is a sum of such terms, the sum of the simplified terms
  1363. is returned.
  1364. Keyword argument ``deep`` controls whether or not subexpressions
  1365. nested deeper inside the main expression are simplified. See examples
  1366. below. Setting `deep` to `False` can save time on nested expressions
  1367. that do not need simplifying on all levels.
  1368. Examples
  1369. ========
  1370. >>> from sympy import symbols
  1371. >>> from sympy.simplify.simplify import nc_simplify
  1372. >>> a, b, c = symbols("a b c", commutative=False)
  1373. >>> nc_simplify(a*b*a*b*c*a*b*c)
  1374. a*b*(a*b*c)**2
  1375. >>> expr = a**2*b*a**4*b*a**4
  1376. >>> nc_simplify(expr)
  1377. a**2*(b*a**4)**2
  1378. >>> nc_simplify(a*b*a*b*c**2*(a*b)**2*c**2)
  1379. ((a*b)**2*c**2)**2
  1380. >>> nc_simplify(a*b*a*b + 2*a*c*a**2*c*a**2*c*a)
  1381. (a*b)**2 + 2*(a*c*a)**3
  1382. >>> nc_simplify(b**-1*a**-1*(a*b)**2)
  1383. a*b
  1384. >>> nc_simplify(a**-1*b**-1*c*a)
  1385. (b*a)**(-1)*c*a
  1386. >>> expr = (a*b*a*b)**2*a*c*a*c
  1387. >>> nc_simplify(expr)
  1388. (a*b)**4*(a*c)**2
  1389. >>> nc_simplify(expr, deep=False)
  1390. (a*b*a*b)**2*(a*c)**2
  1391. '''
  1392. if isinstance(expr, MatrixExpr):
  1393. expr = expr.doit(inv_expand=False)
  1394. _Add, _Mul, _Pow, _Symbol = MatAdd, MatMul, MatPow, MatrixSymbol
  1395. else:
  1396. _Add, _Mul, _Pow, _Symbol = Add, Mul, Pow, Symbol
  1397. # =========== Auxiliary functions ========================
  1398. def _overlaps(args):
  1399. # Calculate a list of lists m such that m[i][j] contains the lengths
  1400. # of all possible overlaps between args[:i+1] and args[i+1+j:].
  1401. # An overlap is a suffix of the prefix that matches a prefix
  1402. # of the suffix.
  1403. # For example, let expr=c*a*b*a*b*a*b*a*b. Then m[3][0] contains
  1404. # the lengths of overlaps of c*a*b*a*b with a*b*a*b. The overlaps
  1405. # are a*b*a*b, a*b and the empty word so that m[3][0]=[4,2,0].
  1406. # All overlaps rather than only the longest one are recorded
  1407. # because this information helps calculate other overlap lengths.
  1408. m = [[([1, 0] if a == args[0] else [0]) for a in args[1:]]]
  1409. for i in range(1, len(args)):
  1410. overlaps = []
  1411. j = 0
  1412. for j in range(len(args) - i - 1):
  1413. overlap = []
  1414. for v in m[i-1][j+1]:
  1415. if j + i + 1 + v < len(args) and args[i] == args[j+i+1+v]:
  1416. overlap.append(v + 1)
  1417. overlap += [0]
  1418. overlaps.append(overlap)
  1419. m.append(overlaps)
  1420. return m
  1421. def _reduce_inverses(_args):
  1422. # replace consecutive negative powers by an inverse
  1423. # of a product of positive powers, e.g. a**-1*b**-1*c
  1424. # will simplify to (a*b)**-1*c;
  1425. # return that new args list and the number of negative
  1426. # powers in it (inv_tot)
  1427. inv_tot = 0 # total number of inverses
  1428. inverses = []
  1429. args = []
  1430. for arg in _args:
  1431. if isinstance(arg, _Pow) and arg.args[1].is_extended_negative:
  1432. inverses = [arg**-1] + inverses
  1433. inv_tot += 1
  1434. else:
  1435. if len(inverses) == 1:
  1436. args.append(inverses[0]**-1)
  1437. elif len(inverses) > 1:
  1438. args.append(_Pow(_Mul(*inverses), -1))
  1439. inv_tot -= len(inverses) - 1
  1440. inverses = []
  1441. args.append(arg)
  1442. if inverses:
  1443. args.append(_Pow(_Mul(*inverses), -1))
  1444. inv_tot -= len(inverses) - 1
  1445. return inv_tot, tuple(args)
  1446. def get_score(s):
  1447. # compute the number of arguments of s
  1448. # (including in nested expressions) overall
  1449. # but ignore exponents
  1450. if isinstance(s, _Pow):
  1451. return get_score(s.args[0])
  1452. elif isinstance(s, (_Add, _Mul)):
  1453. return sum(get_score(a) for a in s.args)
  1454. return 1
  1455. def compare(s, alt_s):
  1456. # compare two possible simplifications and return a
  1457. # "better" one
  1458. if s != alt_s and get_score(alt_s) < get_score(s):
  1459. return alt_s
  1460. return s
  1461. # ========================================================
  1462. if not isinstance(expr, (_Add, _Mul, _Pow)) or expr.is_commutative:
  1463. return expr
  1464. args = expr.args[:]
  1465. if isinstance(expr, _Pow):
  1466. if deep:
  1467. return _Pow(nc_simplify(args[0]), args[1]).doit()
  1468. else:
  1469. return expr
  1470. elif isinstance(expr, _Add):
  1471. return _Add(*[nc_simplify(a, deep=deep) for a in args]).doit()
  1472. else:
  1473. # get the non-commutative part
  1474. c_args, args = expr.args_cnc()
  1475. com_coeff = Mul(*c_args)
  1476. if not equal_valued(com_coeff, 1):
  1477. return com_coeff*nc_simplify(expr/com_coeff, deep=deep)
  1478. inv_tot, args = _reduce_inverses(args)
  1479. # if most arguments are negative, work with the inverse
  1480. # of the expression, e.g. a**-1*b*a**-1*c**-1 will become
  1481. # (c*a*b**-1*a)**-1 at the end so can work with c*a*b**-1*a
  1482. invert = False
  1483. if inv_tot > len(args)/2:
  1484. invert = True
  1485. args = [a**-1 for a in args[::-1]]
  1486. if deep:
  1487. args = tuple(nc_simplify(a) for a in args)
  1488. m = _overlaps(args)
  1489. # simps will be {subterm: end} where `end` is the ending
  1490. # index of a sequence of repetitions of subterm;
  1491. # this is for not wasting time with subterms that are part
  1492. # of longer, already considered sequences
  1493. simps = {}
  1494. post = 1
  1495. pre = 1
  1496. # the simplification coefficient is the number of
  1497. # arguments by which contracting a given sequence
  1498. # would reduce the word; e.g. in a*b*a*b*c*a*b*c,
  1499. # contracting a*b*a*b to (a*b)**2 removes 3 arguments
  1500. # while a*b*c*a*b*c to (a*b*c)**2 removes 6. It's
  1501. # better to contract the latter so simplification
  1502. # with a maximum simplification coefficient will be chosen
  1503. max_simp_coeff = 0
  1504. simp = None # information about future simplification
  1505. for i in range(1, len(args)):
  1506. simp_coeff = 0
  1507. l = 0 # length of a subterm
  1508. p = 0 # the power of a subterm
  1509. if i < len(args) - 1:
  1510. rep = m[i][0]
  1511. start = i # starting index of the repeated sequence
  1512. end = i+1 # ending index of the repeated sequence
  1513. if i == len(args)-1 or rep == [0]:
  1514. # no subterm is repeated at this stage, at least as
  1515. # far as the arguments are concerned - there may be
  1516. # a repetition if powers are taken into account
  1517. if (isinstance(args[i], _Pow) and
  1518. not isinstance(args[i].args[0], _Symbol)):
  1519. subterm = args[i].args[0].args
  1520. l = len(subterm)
  1521. if args[i-l:i] == subterm:
  1522. # e.g. a*b in a*b*(a*b)**2 is not repeated
  1523. # in args (= [a, b, (a*b)**2]) but it
  1524. # can be matched here
  1525. p += 1
  1526. start -= l
  1527. if args[i+1:i+1+l] == subterm:
  1528. # e.g. a*b in (a*b)**2*a*b
  1529. p += 1
  1530. end += l
  1531. if p:
  1532. p += args[i].args[1]
  1533. else:
  1534. continue
  1535. else:
  1536. l = rep[0] # length of the longest repeated subterm at this point
  1537. start -= l - 1
  1538. subterm = args[start:end]
  1539. p = 2
  1540. end += l
  1541. if subterm in simps and simps[subterm] >= start:
  1542. # the subterm is part of a sequence that
  1543. # has already been considered
  1544. continue
  1545. # count how many times it's repeated
  1546. while end < len(args):
  1547. if l in m[end-1][0]:
  1548. p += 1
  1549. end += l
  1550. elif isinstance(args[end], _Pow) and args[end].args[0].args == subterm:
  1551. # for cases like a*b*a*b*(a*b)**2*a*b
  1552. p += args[end].args[1]
  1553. end += 1
  1554. else:
  1555. break
  1556. # see if another match can be made, e.g.
  1557. # for b*a**2 in b*a**2*b*a**3 or a*b in
  1558. # a**2*b*a*b
  1559. pre_exp = 0
  1560. pre_arg = 1
  1561. if start - l >= 0 and args[start-l+1:start] == subterm[1:]:
  1562. if isinstance(subterm[0], _Pow):
  1563. pre_arg = subterm[0].args[0]
  1564. exp = subterm[0].args[1]
  1565. else:
  1566. pre_arg = subterm[0]
  1567. exp = 1
  1568. if isinstance(args[start-l], _Pow) and args[start-l].args[0] == pre_arg:
  1569. pre_exp = args[start-l].args[1] - exp
  1570. start -= l
  1571. p += 1
  1572. elif args[start-l] == pre_arg:
  1573. pre_exp = 1 - exp
  1574. start -= l
  1575. p += 1
  1576. post_exp = 0
  1577. post_arg = 1
  1578. if end + l - 1 < len(args) and args[end:end+l-1] == subterm[:-1]:
  1579. if isinstance(subterm[-1], _Pow):
  1580. post_arg = subterm[-1].args[0]
  1581. exp = subterm[-1].args[1]
  1582. else:
  1583. post_arg = subterm[-1]
  1584. exp = 1
  1585. if isinstance(args[end+l-1], _Pow) and args[end+l-1].args[0] == post_arg:
  1586. post_exp = args[end+l-1].args[1] - exp
  1587. end += l
  1588. p += 1
  1589. elif args[end+l-1] == post_arg:
  1590. post_exp = 1 - exp
  1591. end += l
  1592. p += 1
  1593. # Consider a*b*a**2*b*a**2*b*a:
  1594. # b*a**2 is explicitly repeated, but note
  1595. # that in this case a*b*a is also repeated
  1596. # so there are two possible simplifications:
  1597. # a*(b*a**2)**3*a**-1 or (a*b*a)**3
  1598. # The latter is obviously simpler.
  1599. # But in a*b*a**2*b**2*a**2 the simplifications are
  1600. # a*(b*a**2)**2 and (a*b*a)**3*a in which case
  1601. # it's better to stick with the shorter subterm
  1602. if post_exp and exp % 2 == 0 and start > 0:
  1603. exp = exp/2
  1604. _pre_exp = 1
  1605. _post_exp = 1
  1606. if isinstance(args[start-1], _Pow) and args[start-1].args[0] == post_arg:
  1607. _post_exp = post_exp + exp
  1608. _pre_exp = args[start-1].args[1] - exp
  1609. elif args[start-1] == post_arg:
  1610. _post_exp = post_exp + exp
  1611. _pre_exp = 1 - exp
  1612. if _pre_exp == 0 or _post_exp == 0:
  1613. if not pre_exp:
  1614. start -= 1
  1615. post_exp = _post_exp
  1616. pre_exp = _pre_exp
  1617. pre_arg = post_arg
  1618. subterm = (post_arg**exp,) + subterm[:-1] + (post_arg**exp,)
  1619. simp_coeff += end-start
  1620. if post_exp:
  1621. simp_coeff -= 1
  1622. if pre_exp:
  1623. simp_coeff -= 1
  1624. simps[subterm] = end
  1625. if simp_coeff > max_simp_coeff:
  1626. max_simp_coeff = simp_coeff
  1627. simp = (start, _Mul(*subterm), p, end, l)
  1628. pre = pre_arg**pre_exp
  1629. post = post_arg**post_exp
  1630. if simp:
  1631. subterm = _Pow(nc_simplify(simp[1], deep=deep), simp[2])
  1632. pre = nc_simplify(_Mul(*args[:simp[0]])*pre, deep=deep)
  1633. post = post*nc_simplify(_Mul(*args[simp[3]:]), deep=deep)
  1634. simp = pre*subterm*post
  1635. if pre != 1 or post != 1:
  1636. # new simplifications may be possible but no need
  1637. # to recurse over arguments
  1638. simp = nc_simplify(simp, deep=False)
  1639. else:
  1640. simp = _Mul(*args)
  1641. if invert:
  1642. simp = _Pow(simp, -1)
  1643. # see if factor_nc(expr) is simplified better
  1644. if not isinstance(expr, MatrixExpr):
  1645. f_expr = factor_nc(expr)
  1646. if f_expr != expr:
  1647. alt_simp = nc_simplify(f_expr, deep=deep)
  1648. simp = compare(simp, alt_simp)
  1649. else:
  1650. simp = simp.doit(inv_expand=False)
  1651. return simp
  1652. def dotprodsimp(expr, withsimp=False):
  1653. """Simplification for a sum of products targeted at the kind of blowup that
  1654. occurs during summation of products. Intended to reduce expression blowup
  1655. during matrix multiplication or other similar operations. Only works with
  1656. algebraic expressions and does not recurse into non.
  1657. Parameters
  1658. ==========
  1659. withsimp : bool, optional
  1660. Specifies whether a flag should be returned along with the expression
  1661. to indicate roughly whether simplification was successful. It is used
  1662. in ``MatrixArithmetic._eval_pow_by_recursion`` to avoid attempting to
  1663. simplify an expression repetitively which does not simplify.
  1664. """
  1665. def count_ops_alg(expr):
  1666. """Optimized count algebraic operations with no recursion into
  1667. non-algebraic args that ``core.function.count_ops`` does. Also returns
  1668. whether rational functions may be present according to negative
  1669. exponents of powers or non-number fractions.
  1670. Returns
  1671. =======
  1672. ops, ratfunc : int, bool
  1673. ``ops`` is the number of algebraic operations starting at the top
  1674. level expression (not recursing into non-alg children). ``ratfunc``
  1675. specifies whether the expression MAY contain rational functions
  1676. which ``cancel`` MIGHT optimize.
  1677. """
  1678. ops = 0
  1679. args = [expr]
  1680. ratfunc = False
  1681. while args:
  1682. a = args.pop()
  1683. if not isinstance(a, Basic):
  1684. continue
  1685. if a.is_Rational:
  1686. if a is not S.One: # -1/3 = NEG + DIV
  1687. ops += bool (a.p < 0) + bool (a.q != 1)
  1688. elif a.is_Mul:
  1689. if a.could_extract_minus_sign():
  1690. ops += 1
  1691. if a.args[0] is S.NegativeOne:
  1692. a = a.as_two_terms()[1]
  1693. else:
  1694. a = -a
  1695. n, d = fraction(a)
  1696. if n.is_Integer:
  1697. ops += 1 + bool (n < 0)
  1698. args.append(d) # won't be -Mul but could be Add
  1699. elif d is not S.One:
  1700. if not d.is_Integer:
  1701. args.append(d)
  1702. ratfunc=True
  1703. ops += 1
  1704. args.append(n) # could be -Mul
  1705. else:
  1706. ops += len(a.args) - 1
  1707. args.extend(a.args)
  1708. elif a.is_Add:
  1709. laargs = len(a.args)
  1710. negs = 0
  1711. for ai in a.args:
  1712. if ai.could_extract_minus_sign():
  1713. negs += 1
  1714. ai = -ai
  1715. args.append(ai)
  1716. ops += laargs - (negs != laargs) # -x - y = NEG + SUB
  1717. elif a.is_Pow:
  1718. ops += 1
  1719. args.append(a.base)
  1720. if not ratfunc:
  1721. ratfunc = a.exp.is_negative is not False
  1722. return ops, ratfunc
  1723. def nonalg_subs_dummies(expr, dummies):
  1724. """Substitute dummy variables for non-algebraic expressions to avoid
  1725. evaluation of non-algebraic terms that ``polys.polytools.cancel`` does.
  1726. """
  1727. if not expr.args:
  1728. return expr
  1729. if expr.is_Add or expr.is_Mul or expr.is_Pow:
  1730. args = None
  1731. for i, a in enumerate(expr.args):
  1732. c = nonalg_subs_dummies(a, dummies)
  1733. if c is a:
  1734. continue
  1735. if args is None:
  1736. args = list(expr.args)
  1737. args[i] = c
  1738. if args is None:
  1739. return expr
  1740. return expr.func(*args)
  1741. return dummies.setdefault(expr, Dummy())
  1742. simplified = False # doesn't really mean simplified, rather "can simplify again"
  1743. if isinstance(expr, Basic) and (expr.is_Add or expr.is_Mul or expr.is_Pow):
  1744. expr2 = expr.expand(deep=True, modulus=None, power_base=False,
  1745. power_exp=False, mul=True, log=False, multinomial=True, basic=False)
  1746. if expr2 != expr:
  1747. expr = expr2
  1748. simplified = True
  1749. exprops, ratfunc = count_ops_alg(expr)
  1750. if exprops >= 6: # empirically tested cutoff for expensive simplification
  1751. if ratfunc:
  1752. dummies = {}
  1753. expr2 = nonalg_subs_dummies(expr, dummies)
  1754. if expr2 is expr or count_ops_alg(expr2)[0] >= 6: # check again after substitution
  1755. expr3 = cancel(expr2)
  1756. if expr3 != expr2:
  1757. expr = expr3.subs([(d, e) for e, d in dummies.items()])
  1758. simplified = True
  1759. # very special case: x/(x-1) - 1/(x-1) -> 1
  1760. elif (exprops == 5 and expr.is_Add and expr.args [0].is_Mul and
  1761. expr.args [1].is_Mul and expr.args [0].args [-1].is_Pow and
  1762. expr.args [1].args [-1].is_Pow and
  1763. expr.args [0].args [-1].exp is S.NegativeOne and
  1764. expr.args [1].args [-1].exp is S.NegativeOne):
  1765. expr2 = together (expr)
  1766. expr2ops = count_ops_alg(expr2)[0]
  1767. if expr2ops < exprops:
  1768. expr = expr2
  1769. simplified = True
  1770. else:
  1771. simplified = True
  1772. return (expr, simplified) if withsimp else expr
  1773. bottom_up = deprecated(
  1774. """
  1775. Using bottom_up from the sympy.simplify.simplify submodule is
  1776. deprecated.
  1777. Instead, use bottom_up from the top-level sympy namespace, like
  1778. sympy.bottom_up
  1779. """,
  1780. deprecated_since_version="1.10",
  1781. active_deprecations_target="deprecated-traversal-functions-moved",
  1782. )(_bottom_up)
  1783. # XXX: This function really should either be private API or exported in the
  1784. # top-level sympy/__init__.py
  1785. walk = deprecated(
  1786. """
  1787. Using walk from the sympy.simplify.simplify submodule is
  1788. deprecated.
  1789. Instead, use walk from sympy.core.traversal.walk
  1790. """,
  1791. deprecated_since_version="1.10",
  1792. active_deprecations_target="deprecated-traversal-functions-moved",
  1793. )(_walk)