simplify.py 71 KB

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