Coverage for /usr/lib/python3/dist-packages/sympy/simplify/simplify.py: 14%

977 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1from collections import defaultdict 

2 

3from sympy.concrete.products import Product 

4from sympy.concrete.summations import Sum 

5from sympy.core import (Basic, S, Add, Mul, Pow, Symbol, sympify, 

6 expand_func, Function, Dummy, Expr, factor_terms, 

7 expand_power_exp, Eq) 

8from sympy.core.exprtools import factor_nc 

9from sympy.core.parameters import global_parameters 

10from sympy.core.function import (expand_log, count_ops, _mexpand, 

11 nfloat, expand_mul, expand) 

12from sympy.core.numbers import Float, I, pi, Rational 

13from sympy.core.relational import Relational 

14from sympy.core.rules import Transform 

15from sympy.core.sorting import ordered 

16from sympy.core.sympify import _sympify 

17from sympy.core.traversal import bottom_up as _bottom_up, walk as _walk 

18from sympy.functions import gamma, exp, sqrt, log, exp_polar, re 

19from sympy.functions.combinatorial.factorials import CombinatorialFunction 

20from sympy.functions.elementary.complexes import unpolarify, Abs, sign 

21from sympy.functions.elementary.exponential import ExpBase 

22from sympy.functions.elementary.hyperbolic import HyperbolicFunction 

23from sympy.functions.elementary.integers import ceiling 

24from sympy.functions.elementary.piecewise import (Piecewise, piecewise_fold, 

25 piecewise_simplify) 

26from sympy.functions.elementary.trigonometric import TrigonometricFunction 

27from sympy.functions.special.bessel import (BesselBase, besselj, besseli, 

28 besselk, bessely, jn) 

29from sympy.functions.special.tensor_functions import KroneckerDelta 

30from sympy.integrals.integrals import Integral 

31from sympy.matrices.expressions import (MatrixExpr, MatAdd, MatMul, 

32 MatPow, MatrixSymbol) 

33from sympy.polys import together, cancel, factor 

34from sympy.polys.numberfields.minpoly import _is_sum_surds, _minimal_polynomial_sq 

35from sympy.simplify.combsimp import combsimp 

36from sympy.simplify.cse_opts import sub_pre, sub_post 

37from sympy.simplify.hyperexpand import hyperexpand 

38from sympy.simplify.powsimp import powsimp 

39from sympy.simplify.radsimp import radsimp, fraction, collect_abs 

40from sympy.simplify.sqrtdenest import sqrtdenest 

41from sympy.simplify.trigsimp import trigsimp, exptrigsimp 

42from sympy.utilities.decorator import deprecated 

43from sympy.utilities.iterables import has_variety, sift, subsets, iterable 

44from sympy.utilities.misc import as_int 

45 

46import mpmath 

47 

48 

49def 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 

56 Explanation 

57 =========== 

58 

59 If ``dict=True`` then the separated terms will be returned 

60 in a dictionary keyed to their corresponding symbols. 

61 By default, all symbols in the expression will appear as 

62 keys; if symbols are provided, then all those symbols will 

63 be used as keys, and any terms in the expression containing 

64 other symbols or non-symbols will be returned keyed to the 

65 string 'coeff'. (Passing None for symbols will return the 

66 expression in a dictionary keyed to 'coeff'.) 

67 

68 If ``force=True``, then bases of powers will be separated regardless 

69 of assumptions on the symbols involved. 

70 

71 Notes 

72 ===== 

73 

74 The order of the factors is determined by Mul, so that the 

75 separated expressions may not necessarily be grouped together. 

76 

77 Although factoring is necessary to separate variables in some 

78 expressions, it is not necessary in all cases, so one should not 

79 count on the returned factors being factored. 

80 

81 Examples 

82 ======== 

83 

84 >>> from sympy.abc import x, y, z, alpha 

85 >>> from sympy import separatevars, sin 

86 >>> separatevars((x*y)**y) 

87 (x*y)**y 

88 >>> separatevars((x*y)**y, force=True) 

89 x**y*y**y 

90 

91 >>> e = 2*x**2*z*sin(y)+2*z*x**2 

92 >>> separatevars(e) 

93 2*x**2*z*(sin(y) + 1) 

94 >>> separatevars(e, symbols=(x, y), dict=True) 

95 {'coeff': 2*z, x: x**2, y: sin(y) + 1} 

96 >>> separatevars(e, [x, y, alpha], dict=True) 

97 {'coeff': 2*z, alpha: 1, x: x**2, y: sin(y) + 1} 

98 

99 If the expression is not really separable, or is only partially 

100 separable, separatevars will do the best it can to separate it 

101 by using factoring. 

102 

103 >>> separatevars(x + x*y - 3*x**2) 

104 -x*(3*x - y - 1) 

105 

106 If the expression is not separable then expr is returned unchanged 

107 or (if dict=True) then None is returned. 

108 

109 >>> eq = 2*x + y*sin(x) 

110 >>> separatevars(eq) == eq 

111 True 

112 >>> separatevars(2*x + y*sin(x), symbols=(x, y), dict=True) is None 

113 True 

114 

115 """ 

116 expr = sympify(expr) 

117 if dict: 

118 return _separatevars_dict(_separatevars(expr, force), symbols) 

119 else: 

120 return _separatevars(expr, force) 

121 

122 

123def _separatevars(expr, force): 

124 if isinstance(expr, Abs): 

125 arg = expr.args[0] 

126 if arg.is_Mul and not arg.is_number: 

127 s = separatevars(arg, dict=True, force=force) 

128 if s is not None: 

129 return Mul(*map(expr.func, s.values())) 

130 else: 

131 return expr 

132 

133 if len(expr.free_symbols) < 2: 

134 return expr 

135 

136 # don't destroy a Mul since much of the work may already be done 

137 if expr.is_Mul: 

138 args = list(expr.args) 

139 changed = False 

140 for i, a in enumerate(args): 

141 args[i] = separatevars(a, force) 

142 changed = changed or args[i] != a 

143 if changed: 

144 expr = expr.func(*args) 

145 return expr 

146 

147 # get a Pow ready for expansion 

148 if expr.is_Pow and expr.base != S.Exp1: 

149 expr = Pow(separatevars(expr.base, force=force), expr.exp) 

150 

151 # First try other expansion methods 

152 expr = expr.expand(mul=False, multinomial=False, force=force) 

153 

154 _expr, reps = posify(expr) if force else (expr, {}) 

155 expr = factor(_expr).subs(reps) 

156 

157 if not expr.is_Add: 

158 return expr 

159 

160 # Find any common coefficients to pull out 

161 args = list(expr.args) 

162 commonc = args[0].args_cnc(cset=True, warn=False)[0] 

163 for i in args[1:]: 

164 commonc &= i.args_cnc(cset=True, warn=False)[0] 

165 commonc = Mul(*commonc) 

166 commonc = commonc.as_coeff_Mul()[1] # ignore constants 

167 commonc_set = commonc.args_cnc(cset=True, warn=False)[0] 

168 

169 # remove them 

170 for i, a in enumerate(args): 

171 c, nc = a.args_cnc(cset=True, warn=False) 

172 c = c - commonc_set 

173 args[i] = Mul(*c)*Mul(*nc) 

174 nonsepar = Add(*args) 

175 

176 if len(nonsepar.free_symbols) > 1: 

177 _expr = nonsepar 

178 _expr, reps = posify(_expr) if force else (_expr, {}) 

179 _expr = (factor(_expr)).subs(reps) 

180 

181 if not _expr.is_Add: 

182 nonsepar = _expr 

183 

184 return commonc*nonsepar 

185 

186 

187def _separatevars_dict(expr, symbols): 

188 if symbols: 

189 if not all(t.is_Atom for t in symbols): 

190 raise ValueError("symbols must be Atoms.") 

191 symbols = list(symbols) 

192 elif symbols is None: 

193 return {'coeff': expr} 

194 else: 

195 symbols = list(expr.free_symbols) 

196 if not symbols: 

197 return None 

198 

199 ret = {i: [] for i in symbols + ['coeff']} 

200 

201 for i in Mul.make_args(expr): 

202 expsym = i.free_symbols 

203 intersection = set(symbols).intersection(expsym) 

204 if len(intersection) > 1: 

205 return None 

206 if len(intersection) == 0: 

207 # There are no symbols, so it is part of the coefficient 

208 ret['coeff'].append(i) 

209 else: 

210 ret[intersection.pop()].append(i) 

211 

212 # rebuild 

213 for k, v in ret.items(): 

214 ret[k] = Mul(*v) 

215 

216 return ret 

217 

218 

219def posify(eq): 

220 """Return ``eq`` (with generic symbols made positive) and a 

221 dictionary containing the mapping between the old and new 

222 symbols. 

223 

224 Explanation 

225 =========== 

226 

227 Any symbol that has positive=None will be replaced with a positive dummy 

228 symbol having the same name. This replacement will allow more symbolic 

229 processing of expressions, especially those involving powers and 

230 logarithms. 

231 

232 A dictionary that can be sent to subs to restore ``eq`` to its original 

233 symbols is also returned. 

234 

235 >>> from sympy import posify, Symbol, log, solve 

236 >>> from sympy.abc import x 

237 >>> posify(x + Symbol('p', positive=True) + Symbol('n', negative=True)) 

238 (_x + n + p, {_x: x}) 

239 

240 >>> eq = 1/x 

241 >>> log(eq).expand() 

242 log(1/x) 

243 >>> log(posify(eq)[0]).expand() 

244 -log(_x) 

245 >>> p, rep = posify(eq) 

246 >>> log(p).expand().subs(rep) 

247 -log(x) 

248 

249 It is possible to apply the same transformations to an iterable 

250 of expressions: 

251 

252 >>> eq = x**2 - 4 

253 >>> solve(eq, x) 

254 [-2, 2] 

255 >>> eq_x, reps = posify([eq, x]); eq_x 

256 [_x**2 - 4, _x] 

257 >>> solve(*eq_x) 

258 [2] 

259 """ 

260 eq = sympify(eq) 

261 if iterable(eq): 

262 f = type(eq) 

263 eq = list(eq) 

264 syms = set() 

265 for e in eq: 

266 syms = syms.union(e.atoms(Symbol)) 

267 reps = {} 

268 for s in syms: 

269 reps.update({v: k for k, v in posify(s)[1].items()}) 

270 for i, e in enumerate(eq): 

271 eq[i] = e.subs(reps) 

272 return f(eq), {r: s for s, r in reps.items()} 

273 

274 reps = {s: Dummy(s.name, positive=True, **s.assumptions0) 

275 for s in eq.free_symbols if s.is_positive is None} 

276 eq = eq.subs(reps) 

277 return eq, {r: s for s, r in reps.items()} 

278 

279 

280def hypersimp(f, k): 

281 """Given combinatorial term f(k) simplify its consecutive term ratio 

282 i.e. f(k+1)/f(k). The input term can be composed of functions and 

283 integer sequences which have equivalent representation in terms 

284 of gamma special function. 

285 

286 Explanation 

287 =========== 

288 

289 The algorithm performs three basic steps: 

290 

291 1. Rewrite all functions in terms of gamma, if possible. 

292 

293 2. Rewrite all occurrences of gamma in terms of products 

294 of gamma and rising factorial with integer, absolute 

295 constant exponent. 

296 

297 3. Perform simplification of nested fractions, powers 

298 and if the resulting expression is a quotient of 

299 polynomials, reduce their total degree. 

300 

301 If f(k) is hypergeometric then as result we arrive with a 

302 quotient of polynomials of minimal degree. Otherwise None 

303 is returned. 

304 

305 For more information on the implemented algorithm refer to: 

306 

307 1. W. Koepf, Algorithms for m-fold Hypergeometric Summation, 

308 Journal of Symbolic Computation (1995) 20, 399-417 

309 """ 

310 f = sympify(f) 

311 

312 g = f.subs(k, k + 1) / f 

313 

314 g = g.rewrite(gamma) 

315 if g.has(Piecewise): 

316 g = piecewise_fold(g) 

317 g = g.args[-1][0] 

318 g = expand_func(g) 

319 g = powsimp(g, deep=True, combine='exp') 

320 

321 if g.is_rational_function(k): 

322 return simplify(g, ratio=S.Infinity) 

323 else: 

324 return None 

325 

326 

327def hypersimilar(f, g, k): 

328 """ 

329 Returns True if ``f`` and ``g`` are hyper-similar. 

330 

331 Explanation 

332 =========== 

333 

334 Similarity in hypergeometric sense means that a quotient of 

335 f(k) and g(k) is a rational function in ``k``. This procedure 

336 is useful in solving recurrence relations. 

337 

338 For more information see hypersimp(). 

339 

340 """ 

341 f, g = list(map(sympify, (f, g))) 

342 

343 h = (f/g).rewrite(gamma) 

344 h = h.expand(func=True, basic=False) 

345 

346 return h.is_rational_function(k) 

347 

348 

349def signsimp(expr, evaluate=None): 

350 """Make all Add sub-expressions canonical wrt sign. 

351 

352 Explanation 

353 =========== 

354 

355 If an Add subexpression, ``a``, can have a sign extracted, 

356 as determined by could_extract_minus_sign, it is replaced 

357 with Mul(-1, a, evaluate=False). This allows signs to be 

358 extracted from powers and products. 

359 

360 Examples 

361 ======== 

362 

363 >>> from sympy import signsimp, exp, symbols 

364 >>> from sympy.abc import x, y 

365 >>> i = symbols('i', odd=True) 

366 >>> n = -1 + 1/x 

367 >>> n/x/(-n)**2 - 1/n/x 

368 (-1 + 1/x)/(x*(1 - 1/x)**2) - 1/(x*(-1 + 1/x)) 

369 >>> signsimp(_) 

370 0 

371 >>> x*n + x*-n 

372 x*(-1 + 1/x) + x*(1 - 1/x) 

373 >>> signsimp(_) 

374 0 

375 

376 Since powers automatically handle leading signs 

377 

378 >>> (-2)**i 

379 -2**i 

380 

381 signsimp can be used to put the base of a power with an integer 

382 exponent into canonical form: 

383 

384 >>> n**i 

385 (-1 + 1/x)**i 

386 

387 By default, signsimp does not leave behind any hollow simplification: 

388 if making an Add canonical wrt sign didn't change the expression, the 

389 original Add is restored. If this is not desired then the keyword 

390 ``evaluate`` can be set to False: 

391 

392 >>> e = exp(y - x) 

393 >>> signsimp(e) == e 

394 True 

395 >>> signsimp(e, evaluate=False) 

396 exp(-(x - y)) 

397 

398 """ 

399 if evaluate is None: 

400 evaluate = global_parameters.evaluate 

401 expr = sympify(expr) 

402 if not isinstance(expr, (Expr, Relational)) or expr.is_Atom: 

403 return expr 

404 # get rid of an pre-existing unevaluation regarding sign 

405 e = expr.replace(lambda x: x.is_Mul and -(-x) != x, lambda x: -(-x)) 

406 e = sub_post(sub_pre(e)) 

407 if not isinstance(e, (Expr, Relational)) or e.is_Atom: 

408 return e 

409 if e.is_Add: 

410 rv = e.func(*[signsimp(a) for a in e.args]) 

411 if not evaluate and isinstance(rv, Add 

412 ) and rv.could_extract_minus_sign(): 

413 return Mul(S.NegativeOne, -rv, evaluate=False) 

414 return rv 

415 if evaluate: 

416 e = e.replace(lambda x: x.is_Mul and -(-x) != x, lambda x: -(-x)) 

417 return e 

418 

419 

420def simplify(expr, ratio=1.7, measure=count_ops, rational=False, inverse=False, doit=True, **kwargs): 

421 """Simplifies the given expression. 

422 

423 Explanation 

424 =========== 

425 

426 Simplification is not a well defined term and the exact strategies 

427 this function tries can change in the future versions of SymPy. If 

428 your algorithm relies on "simplification" (whatever it is), try to 

429 determine what you need exactly - is it powsimp()?, radsimp()?, 

430 together()?, logcombine()?, or something else? And use this particular 

431 function directly, because those are well defined and thus your algorithm 

432 will be robust. 

433 

434 Nonetheless, especially for interactive use, or when you do not know 

435 anything about the structure of the expression, simplify() tries to apply 

436 intelligent heuristics to make the input expression "simpler". For 

437 example: 

438 

439 >>> from sympy import simplify, cos, sin 

440 >>> from sympy.abc import x, y 

441 >>> a = (x + x**2)/(x*sin(y)**2 + x*cos(y)**2) 

442 >>> a 

443 (x**2 + x)/(x*sin(y)**2 + x*cos(y)**2) 

444 >>> simplify(a) 

445 x + 1 

446 

447 Note that we could have obtained the same result by using specific 

448 simplification functions: 

449 

450 >>> from sympy import trigsimp, cancel 

451 >>> trigsimp(a) 

452 (x**2 + x)/x 

453 >>> cancel(_) 

454 x + 1 

455 

456 In some cases, applying :func:`simplify` may actually result in some more 

457 complicated expression. The default ``ratio=1.7`` prevents more extreme 

458 cases: if (result length)/(input length) > ratio, then input is returned 

459 unmodified. The ``measure`` parameter lets you specify the function used 

460 to determine how complex an expression is. The function should take a 

461 single argument as an expression and return a number such that if 

462 expression ``a`` is more complex than expression ``b``, then 

463 ``measure(a) > measure(b)``. The default measure function is 

464 :func:`~.count_ops`, which returns the total number of operations in the 

465 expression. 

466 

467 For example, if ``ratio=1``, ``simplify`` output cannot be longer 

468 than input. 

469 

470 :: 

471 

472 >>> from sympy import sqrt, simplify, count_ops, oo 

473 >>> root = 1/(sqrt(2)+3) 

474 

475 Since ``simplify(root)`` would result in a slightly longer expression, 

476 root is returned unchanged instead:: 

477 

478 >>> simplify(root, ratio=1) == root 

479 True 

480 

481 If ``ratio=oo``, simplify will be applied anyway:: 

482 

483 >>> count_ops(simplify(root, ratio=oo)) > count_ops(root) 

484 True 

485 

486 Note that the shortest expression is not necessary the simplest, so 

487 setting ``ratio`` to 1 may not be a good idea. 

488 Heuristically, the default value ``ratio=1.7`` seems like a reasonable 

489 choice. 

490 

491 You can easily define your own measure function based on what you feel 

492 should represent the "size" or "complexity" of the input expression. Note 

493 that some choices, such as ``lambda expr: len(str(expr))`` may appear to be 

494 good metrics, but have other problems (in this case, the measure function 

495 may slow down simplify too much for very large expressions). If you do not 

496 know what a good metric would be, the default, ``count_ops``, is a good 

497 one. 

498 

499 For example: 

500 

501 >>> from sympy import symbols, log 

502 >>> a, b = symbols('a b', positive=True) 

503 >>> g = log(a) + log(b) + log(a)*log(1/b) 

504 >>> h = simplify(g) 

505 >>> h 

506 log(a*b**(1 - log(a))) 

507 >>> count_ops(g) 

508 8 

509 >>> count_ops(h) 

510 5 

511 

512 So you can see that ``h`` is simpler than ``g`` using the count_ops metric. 

513 However, we may not like how ``simplify`` (in this case, using 

514 ``logcombine``) has created the ``b**(log(1/a) + 1)`` term. A simple way 

515 to reduce this would be to give more weight to powers as operations in 

516 ``count_ops``. We can do this by using the ``visual=True`` option: 

517 

518 >>> print(count_ops(g, visual=True)) 

519 2*ADD + DIV + 4*LOG + MUL 

520 >>> print(count_ops(h, visual=True)) 

521 2*LOG + MUL + POW + SUB 

522 

523 >>> from sympy import Symbol, S 

524 >>> def my_measure(expr): 

525 ... POW = Symbol('POW') 

526 ... # Discourage powers by giving POW a weight of 10 

527 ... count = count_ops(expr, visual=True).subs(POW, 10) 

528 ... # Every other operation gets a weight of 1 (the default) 

529 ... count = count.replace(Symbol, type(S.One)) 

530 ... return count 

531 >>> my_measure(g) 

532 8 

533 >>> my_measure(h) 

534 14 

535 >>> 15./8 > 1.7 # 1.7 is the default ratio 

536 True 

537 >>> simplify(g, measure=my_measure) 

538 -log(a)*log(b) + log(a) + log(b) 

539 

540 Note that because ``simplify()`` internally tries many different 

541 simplification strategies and then compares them using the measure 

542 function, we get a completely different result that is still different 

543 from the input expression by doing this. 

544 

545 If ``rational=True``, Floats will be recast as Rationals before simplification. 

546 If ``rational=None``, Floats will be recast as Rationals but the result will 

547 be recast as Floats. If rational=False(default) then nothing will be done 

548 to the Floats. 

549 

550 If ``inverse=True``, it will be assumed that a composition of inverse 

551 functions, such as sin and asin, can be cancelled in any order. 

552 For example, ``asin(sin(x))`` will yield ``x`` without checking whether 

553 x belongs to the set where this relation is true. The default is 

554 False. 

555 

556 Note that ``simplify()`` automatically calls ``doit()`` on the final 

557 expression. You can avoid this behavior by passing ``doit=False`` as 

558 an argument. 

559 

560 Also, it should be noted that simplifying a boolean expression is not 

561 well defined. If the expression prefers automatic evaluation (such as 

562 :obj:`~.Eq()` or :obj:`~.Or()`), simplification will return ``True`` or 

563 ``False`` if truth value can be determined. If the expression is not 

564 evaluated by default (such as :obj:`~.Predicate()`), simplification will 

565 not reduce it and you should use :func:`~.refine()` or :func:`~.ask()` 

566 function. This inconsistency will be resolved in future version. 

567 

568 See Also 

569 ======== 

570 

571 sympy.assumptions.refine.refine : Simplification using assumptions. 

572 sympy.assumptions.ask.ask : Query for boolean expressions using assumptions. 

573 """ 

574 

575 def shorter(*choices): 

576 """ 

577 Return the choice that has the fewest ops. In case of a tie, 

578 the expression listed first is selected. 

579 """ 

580 if not has_variety(choices): 

581 return choices[0] 

582 return min(choices, key=measure) 

583 

584 def done(e): 

585 rv = e.doit() if doit else e 

586 return shorter(rv, collect_abs(rv)) 

587 

588 expr = sympify(expr, rational=rational) 

589 kwargs = { 

590 "ratio": kwargs.get('ratio', ratio), 

591 "measure": kwargs.get('measure', measure), 

592 "rational": kwargs.get('rational', rational), 

593 "inverse": kwargs.get('inverse', inverse), 

594 "doit": kwargs.get('doit', doit)} 

595 # no routine for Expr needs to check for is_zero 

596 if isinstance(expr, Expr) and expr.is_zero: 

597 return S.Zero if not expr.is_Number else expr 

598 

599 _eval_simplify = getattr(expr, '_eval_simplify', None) 

600 if _eval_simplify is not None: 

601 return _eval_simplify(**kwargs) 

602 

603 original_expr = expr = collect_abs(signsimp(expr)) 

604 

605 if not isinstance(expr, Basic) or not expr.args: # XXX: temporary hack 

606 return expr 

607 

608 if inverse and expr.has(Function): 

609 expr = inversecombine(expr) 

610 if not expr.args: # simplified to atomic 

611 return expr 

612 

613 # do deep simplification 

614 handled = Add, Mul, Pow, ExpBase 

615 expr = expr.replace( 

616 # here, checking for x.args is not enough because Basic has 

617 # args but Basic does not always play well with replace, e.g. 

618 # when simultaneous is True found expressions will be masked 

619 # off with a Dummy but not all Basic objects in an expression 

620 # can be replaced with a Dummy 

621 lambda x: isinstance(x, Expr) and x.args and not isinstance( 

622 x, handled), 

623 lambda x: x.func(*[simplify(i, **kwargs) for i in x.args]), 

624 simultaneous=False) 

625 if not isinstance(expr, handled): 

626 return done(expr) 

627 

628 if not expr.is_commutative: 

629 expr = nc_simplify(expr) 

630 

631 # TODO: Apply different strategies, considering expression pattern: 

632 # is it a purely rational function? Is there any trigonometric function?... 

633 # See also https://github.com/sympy/sympy/pull/185. 

634 

635 

636 # rationalize Floats 

637 floats = False 

638 if rational is not False and expr.has(Float): 

639 floats = True 

640 expr = nsimplify(expr, rational=True) 

641 

642 expr = _bottom_up(expr, lambda w: getattr(w, 'normal', lambda: w)()) 

643 expr = Mul(*powsimp(expr).as_content_primitive()) 

644 _e = cancel(expr) 

645 expr1 = shorter(_e, _mexpand(_e).cancel()) # issue 6829 

646 expr2 = shorter(together(expr, deep=True), together(expr1, deep=True)) 

647 

648 if ratio is S.Infinity: 

649 expr = expr2 

650 else: 

651 expr = shorter(expr2, expr1, expr) 

652 if not isinstance(expr, Basic): # XXX: temporary hack 

653 return expr 

654 

655 expr = factor_terms(expr, sign=False) 

656 

657 # must come before `Piecewise` since this introduces more `Piecewise` terms 

658 if expr.has(sign): 

659 expr = expr.rewrite(Abs) 

660 

661 # Deal with Piecewise separately to avoid recursive growth of expressions 

662 if expr.has(Piecewise): 

663 # Fold into a single Piecewise 

664 expr = piecewise_fold(expr) 

665 # Apply doit, if doit=True 

666 expr = done(expr) 

667 # Still a Piecewise? 

668 if expr.has(Piecewise): 

669 # Fold into a single Piecewise, in case doit lead to some 

670 # expressions being Piecewise 

671 expr = piecewise_fold(expr) 

672 # kroneckersimp also affects Piecewise 

673 if expr.has(KroneckerDelta): 

674 expr = kroneckersimp(expr) 

675 # Still a Piecewise? 

676 if expr.has(Piecewise): 

677 # Do not apply doit on the segments as it has already 

678 # been done above, but simplify 

679 expr = piecewise_simplify(expr, deep=True, doit=False) 

680 # Still a Piecewise? 

681 if expr.has(Piecewise): 

682 # Try factor common terms 

683 expr = shorter(expr, factor_terms(expr)) 

684 # As all expressions have been simplified above with the 

685 # complete simplify, nothing more needs to be done here 

686 return expr 

687 

688 # hyperexpand automatically only works on hypergeometric terms 

689 # Do this after the Piecewise part to avoid recursive expansion 

690 expr = hyperexpand(expr) 

691 

692 if expr.has(KroneckerDelta): 

693 expr = kroneckersimp(expr) 

694 

695 if expr.has(BesselBase): 

696 expr = besselsimp(expr) 

697 

698 if expr.has(TrigonometricFunction, HyperbolicFunction): 

699 expr = trigsimp(expr, deep=True) 

700 

701 if expr.has(log): 

702 expr = shorter(expand_log(expr, deep=True), logcombine(expr)) 

703 

704 if expr.has(CombinatorialFunction, gamma): 

705 # expression with gamma functions or non-integer arguments is 

706 # automatically passed to gammasimp 

707 expr = combsimp(expr) 

708 

709 if expr.has(Sum): 

710 expr = sum_simplify(expr, **kwargs) 

711 

712 if expr.has(Integral): 

713 expr = expr.xreplace({ 

714 i: factor_terms(i) for i in expr.atoms(Integral)}) 

715 

716 if expr.has(Product): 

717 expr = product_simplify(expr, **kwargs) 

718 

719 from sympy.physics.units import Quantity 

720 

721 if expr.has(Quantity): 

722 from sympy.physics.units.util import quantity_simplify 

723 expr = quantity_simplify(expr) 

724 

725 short = shorter(powsimp(expr, combine='exp', deep=True), powsimp(expr), expr) 

726 short = shorter(short, cancel(short)) 

727 short = shorter(short, factor_terms(short), expand_power_exp(expand_mul(short))) 

728 if short.has(TrigonometricFunction, HyperbolicFunction, ExpBase, exp): 

729 short = exptrigsimp(short) 

730 

731 # get rid of hollow 2-arg Mul factorization 

732 hollow_mul = Transform( 

733 lambda x: Mul(*x.args), 

734 lambda x: 

735 x.is_Mul and 

736 len(x.args) == 2 and 

737 x.args[0].is_Number and 

738 x.args[1].is_Add and 

739 x.is_commutative) 

740 expr = short.xreplace(hollow_mul) 

741 

742 numer, denom = expr.as_numer_denom() 

743 if denom.is_Add: 

744 n, d = fraction(radsimp(1/denom, symbolic=False, max_terms=1)) 

745 if n is not S.One: 

746 expr = (numer*n).expand()/d 

747 

748 if expr.could_extract_minus_sign(): 

749 n, d = fraction(expr) 

750 if d != 0: 

751 expr = signsimp(-n/(-d)) 

752 

753 if measure(expr) > ratio*measure(original_expr): 

754 expr = original_expr 

755 

756 # restore floats 

757 if floats and rational is None: 

758 expr = nfloat(expr, exponent=False) 

759 

760 return done(expr) 

761 

762 

763def sum_simplify(s, **kwargs): 

764 """Main function for Sum simplification""" 

765 if not isinstance(s, Add): 

766 s = s.xreplace({a: sum_simplify(a, **kwargs) 

767 for a in s.atoms(Add) if a.has(Sum)}) 

768 s = expand(s) 

769 if not isinstance(s, Add): 

770 return s 

771 

772 terms = s.args 

773 s_t = [] # Sum Terms 

774 o_t = [] # Other Terms 

775 

776 for term in terms: 

777 sum_terms, other = sift(Mul.make_args(term), 

778 lambda i: isinstance(i, Sum), binary=True) 

779 if not sum_terms: 

780 o_t.append(term) 

781 continue 

782 other = [Mul(*other)] 

783 s_t.append(Mul(*(other + [s._eval_simplify(**kwargs) for s in sum_terms]))) 

784 

785 result = Add(sum_combine(s_t), *o_t) 

786 

787 return result 

788 

789 

790def sum_combine(s_t): 

791 """Helper function for Sum simplification 

792 

793 Attempts to simplify a list of sums, by combining limits / sum function's 

794 returns the simplified sum 

795 """ 

796 used = [False] * len(s_t) 

797 

798 for method in range(2): 

799 for i, s_term1 in enumerate(s_t): 

800 if not used[i]: 

801 for j, s_term2 in enumerate(s_t): 

802 if not used[j] and i != j: 

803 temp = sum_add(s_term1, s_term2, method) 

804 if isinstance(temp, (Sum, Mul)): 

805 s_t[i] = temp 

806 s_term1 = s_t[i] 

807 used[j] = True 

808 

809 result = S.Zero 

810 for i, s_term in enumerate(s_t): 

811 if not used[i]: 

812 result = Add(result, s_term) 

813 

814 return result 

815 

816 

817def factor_sum(self, limits=None, radical=False, clear=False, fraction=False, sign=True): 

818 """Return Sum with constant factors extracted. 

819 

820 If ``limits`` is specified then ``self`` is the summand; the other 

821 keywords are passed to ``factor_terms``. 

822 

823 Examples 

824 ======== 

825 

826 >>> from sympy import Sum 

827 >>> from sympy.abc import x, y 

828 >>> from sympy.simplify.simplify import factor_sum 

829 >>> s = Sum(x*y, (x, 1, 3)) 

830 >>> factor_sum(s) 

831 y*Sum(x, (x, 1, 3)) 

832 >>> factor_sum(s.function, s.limits) 

833 y*Sum(x, (x, 1, 3)) 

834 """ 

835 # XXX deprecate in favor of direct call to factor_terms 

836 kwargs = {"radical": radical, "clear": clear, 

837 "fraction": fraction, "sign": sign} 

838 expr = Sum(self, *limits) if limits else self 

839 return factor_terms(expr, **kwargs) 

840 

841 

842def sum_add(self, other, method=0): 

843 """Helper function for Sum simplification""" 

844 #we know this is something in terms of a constant * a sum 

845 #so we temporarily put the constants inside for simplification 

846 #then simplify the result 

847 def __refactor(val): 

848 args = Mul.make_args(val) 

849 sumv = next(x for x in args if isinstance(x, Sum)) 

850 constant = Mul(*[x for x in args if x != sumv]) 

851 return Sum(constant * sumv.function, *sumv.limits) 

852 

853 if isinstance(self, Mul): 

854 rself = __refactor(self) 

855 else: 

856 rself = self 

857 

858 if isinstance(other, Mul): 

859 rother = __refactor(other) 

860 else: 

861 rother = other 

862 

863 if type(rself) is type(rother): 

864 if method == 0: 

865 if rself.limits == rother.limits: 

866 return factor_sum(Sum(rself.function + rother.function, *rself.limits)) 

867 elif method == 1: 

868 if simplify(rself.function - rother.function) == 0: 

869 if len(rself.limits) == len(rother.limits) == 1: 

870 i = rself.limits[0][0] 

871 x1 = rself.limits[0][1] 

872 y1 = rself.limits[0][2] 

873 j = rother.limits[0][0] 

874 x2 = rother.limits[0][1] 

875 y2 = rother.limits[0][2] 

876 

877 if i == j: 

878 if x2 == y1 + 1: 

879 return factor_sum(Sum(rself.function, (i, x1, y2))) 

880 elif x1 == y2 + 1: 

881 return factor_sum(Sum(rself.function, (i, x2, y1))) 

882 

883 return Add(self, other) 

884 

885 

886def product_simplify(s, **kwargs): 

887 """Main function for Product simplification""" 

888 terms = Mul.make_args(s) 

889 p_t = [] # Product Terms 

890 o_t = [] # Other Terms 

891 

892 deep = kwargs.get('deep', True) 

893 for term in terms: 

894 if isinstance(term, Product): 

895 if deep: 

896 p_t.append(Product(term.function.simplify(**kwargs), 

897 *term.limits)) 

898 else: 

899 p_t.append(term) 

900 else: 

901 o_t.append(term) 

902 

903 used = [False] * len(p_t) 

904 

905 for method in range(2): 

906 for i, p_term1 in enumerate(p_t): 

907 if not used[i]: 

908 for j, p_term2 in enumerate(p_t): 

909 if not used[j] and i != j: 

910 tmp_prod = product_mul(p_term1, p_term2, method) 

911 if isinstance(tmp_prod, Product): 

912 p_t[i] = tmp_prod 

913 used[j] = True 

914 

915 result = Mul(*o_t) 

916 

917 for i, p_term in enumerate(p_t): 

918 if not used[i]: 

919 result = Mul(result, p_term) 

920 

921 return result 

922 

923 

924def product_mul(self, other, method=0): 

925 """Helper function for Product simplification""" 

926 if type(self) is type(other): 

927 if method == 0: 

928 if self.limits == other.limits: 

929 return Product(self.function * other.function, *self.limits) 

930 elif method == 1: 

931 if simplify(self.function - other.function) == 0: 

932 if len(self.limits) == len(other.limits) == 1: 

933 i = self.limits[0][0] 

934 x1 = self.limits[0][1] 

935 y1 = self.limits[0][2] 

936 j = other.limits[0][0] 

937 x2 = other.limits[0][1] 

938 y2 = other.limits[0][2] 

939 

940 if i == j: 

941 if x2 == y1 + 1: 

942 return Product(self.function, (i, x1, y2)) 

943 elif x1 == y2 + 1: 

944 return Product(self.function, (i, x2, y1)) 

945 

946 return Mul(self, other) 

947 

948 

949def _nthroot_solve(p, n, prec): 

950 """ 

951 helper function for ``nthroot`` 

952 It denests ``p**Rational(1, n)`` using its minimal polynomial 

953 """ 

954 from sympy.solvers import solve 

955 while n % 2 == 0: 

956 p = sqrtdenest(sqrt(p)) 

957 n = n // 2 

958 if n == 1: 

959 return p 

960 pn = p**Rational(1, n) 

961 x = Symbol('x') 

962 f = _minimal_polynomial_sq(p, n, x) 

963 if f is None: 

964 return None 

965 sols = solve(f, x) 

966 for sol in sols: 

967 if abs(sol - pn).n() < 1./10**prec: 

968 sol = sqrtdenest(sol) 

969 if _mexpand(sol**n) == p: 

970 return sol 

971 

972 

973def logcombine(expr, force=False): 

974 """ 

975 Takes logarithms and combines them using the following rules: 

976 

977 - log(x) + log(y) == log(x*y) if both are positive 

978 - a*log(x) == log(x**a) if x is positive and a is real 

979 

980 If ``force`` is ``True`` then the assumptions above will be assumed to hold if 

981 there is no assumption already in place on a quantity. For example, if 

982 ``a`` is imaginary or the argument negative, force will not perform a 

983 combination but if ``a`` is a symbol with no assumptions the change will 

984 take place. 

985 

986 Examples 

987 ======== 

988 

989 >>> from sympy import Symbol, symbols, log, logcombine, I 

990 >>> from sympy.abc import a, x, y, z 

991 >>> logcombine(a*log(x) + log(y) - log(z)) 

992 a*log(x) + log(y) - log(z) 

993 >>> logcombine(a*log(x) + log(y) - log(z), force=True) 

994 log(x**a*y/z) 

995 >>> x,y,z = symbols('x,y,z', positive=True) 

996 >>> a = Symbol('a', real=True) 

997 >>> logcombine(a*log(x) + log(y) - log(z)) 

998 log(x**a*y/z) 

999 

1000 The transformation is limited to factors and/or terms that 

1001 contain logs, so the result depends on the initial state of 

1002 expansion: 

1003 

1004 >>> eq = (2 + 3*I)*log(x) 

1005 >>> logcombine(eq, force=True) == eq 

1006 True 

1007 >>> logcombine(eq.expand(), force=True) 

1008 log(x**2) + I*log(x**3) 

1009 

1010 See Also 

1011 ======== 

1012 

1013 posify: replace all symbols with symbols having positive assumptions 

1014 sympy.core.function.expand_log: expand the logarithms of products 

1015 and powers; the opposite of logcombine 

1016 

1017 """ 

1018 

1019 def f(rv): 

1020 if not (rv.is_Add or rv.is_Mul): 

1021 return rv 

1022 

1023 def gooda(a): 

1024 # bool to tell whether the leading ``a`` in ``a*log(x)`` 

1025 # could appear as log(x**a) 

1026 return (a is not S.NegativeOne and # -1 *could* go, but we disallow 

1027 (a.is_extended_real or force and a.is_extended_real is not False)) 

1028 

1029 def goodlog(l): 

1030 # bool to tell whether log ``l``'s argument can combine with others 

1031 a = l.args[0] 

1032 return a.is_positive or force and a.is_nonpositive is not False 

1033 

1034 other = [] 

1035 logs = [] 

1036 log1 = defaultdict(list) 

1037 for a in Add.make_args(rv): 

1038 if isinstance(a, log) and goodlog(a): 

1039 log1[()].append(([], a)) 

1040 elif not a.is_Mul: 

1041 other.append(a) 

1042 else: 

1043 ot = [] 

1044 co = [] 

1045 lo = [] 

1046 for ai in a.args: 

1047 if ai.is_Rational and ai < 0: 

1048 ot.append(S.NegativeOne) 

1049 co.append(-ai) 

1050 elif isinstance(ai, log) and goodlog(ai): 

1051 lo.append(ai) 

1052 elif gooda(ai): 

1053 co.append(ai) 

1054 else: 

1055 ot.append(ai) 

1056 if len(lo) > 1: 

1057 logs.append((ot, co, lo)) 

1058 elif lo: 

1059 log1[tuple(ot)].append((co, lo[0])) 

1060 else: 

1061 other.append(a) 

1062 

1063 # if there is only one log in other, put it with the 

1064 # good logs 

1065 if len(other) == 1 and isinstance(other[0], log): 

1066 log1[()].append(([], other.pop())) 

1067 # if there is only one log at each coefficient and none have 

1068 # an exponent to place inside the log then there is nothing to do 

1069 if not logs and all(len(log1[k]) == 1 and log1[k][0] == [] for k in log1): 

1070 return rv 

1071 

1072 # collapse multi-logs as far as possible in a canonical way 

1073 # TODO: see if x*log(a)+x*log(a)*log(b) -> x*log(a)*(1+log(b))? 

1074 # -- in this case, it's unambiguous, but if it were were a log(c) in 

1075 # each term then it's arbitrary whether they are grouped by log(a) or 

1076 # by log(c). So for now, just leave this alone; it's probably better to 

1077 # let the user decide 

1078 for o, e, l in logs: 

1079 l = list(ordered(l)) 

1080 e = log(l.pop(0).args[0]**Mul(*e)) 

1081 while l: 

1082 li = l.pop(0) 

1083 e = log(li.args[0]**e) 

1084 c, l = Mul(*o), e 

1085 if isinstance(l, log): # it should be, but check to be sure 

1086 log1[(c,)].append(([], l)) 

1087 else: 

1088 other.append(c*l) 

1089 

1090 # logs that have the same coefficient can multiply 

1091 for k in list(log1.keys()): 

1092 log1[Mul(*k)] = log(logcombine(Mul(*[ 

1093 l.args[0]**Mul(*c) for c, l in log1.pop(k)]), 

1094 force=force), evaluate=False) 

1095 

1096 # logs that have oppositely signed coefficients can divide 

1097 for k in ordered(list(log1.keys())): 

1098 if k not in log1: # already popped as -k 

1099 continue 

1100 if -k in log1: 

1101 # figure out which has the minus sign; the one with 

1102 # more op counts should be the one 

1103 num, den = k, -k 

1104 if num.count_ops() > den.count_ops(): 

1105 num, den = den, num 

1106 other.append( 

1107 num*log(log1.pop(num).args[0]/log1.pop(den).args[0], 

1108 evaluate=False)) 

1109 else: 

1110 other.append(k*log1.pop(k)) 

1111 

1112 return Add(*other) 

1113 

1114 return _bottom_up(expr, f) 

1115 

1116 

1117def inversecombine(expr): 

1118 """Simplify the composition of a function and its inverse. 

1119 

1120 Explanation 

1121 =========== 

1122 

1123 No attention is paid to whether the inverse is a left inverse or a 

1124 right inverse; thus, the result will in general not be equivalent 

1125 to the original expression. 

1126 

1127 Examples 

1128 ======== 

1129 

1130 >>> from sympy.simplify.simplify import inversecombine 

1131 >>> from sympy import asin, sin, log, exp 

1132 >>> from sympy.abc import x 

1133 >>> inversecombine(asin(sin(x))) 

1134 x 

1135 >>> inversecombine(2*log(exp(3*x))) 

1136 6*x 

1137 """ 

1138 

1139 def f(rv): 

1140 if isinstance(rv, log): 

1141 if isinstance(rv.args[0], exp) or (rv.args[0].is_Pow and rv.args[0].base == S.Exp1): 

1142 rv = rv.args[0].exp 

1143 elif rv.is_Function and hasattr(rv, "inverse"): 

1144 if (len(rv.args) == 1 and len(rv.args[0].args) == 1 and 

1145 isinstance(rv.args[0], rv.inverse(argindex=1))): 

1146 rv = rv.args[0].args[0] 

1147 if rv.is_Pow and rv.base == S.Exp1: 

1148 if isinstance(rv.exp, log): 

1149 rv = rv.exp.args[0] 

1150 return rv 

1151 

1152 return _bottom_up(expr, f) 

1153 

1154 

1155def kroneckersimp(expr): 

1156 """ 

1157 Simplify expressions with KroneckerDelta. 

1158 

1159 The only simplification currently attempted is to identify multiplicative cancellation: 

1160 

1161 Examples 

1162 ======== 

1163 

1164 >>> from sympy import KroneckerDelta, kroneckersimp 

1165 >>> from sympy.abc import i 

1166 >>> kroneckersimp(1 + KroneckerDelta(0, i) * KroneckerDelta(1, i)) 

1167 1 

1168 """ 

1169 def args_cancel(args1, args2): 

1170 for i1 in range(2): 

1171 for i2 in range(2): 

1172 a1 = args1[i1] 

1173 a2 = args2[i2] 

1174 a3 = args1[(i1 + 1) % 2] 

1175 a4 = args2[(i2 + 1) % 2] 

1176 if Eq(a1, a2) is S.true and Eq(a3, a4) is S.false: 

1177 return True 

1178 return False 

1179 

1180 def cancel_kronecker_mul(m): 

1181 args = m.args 

1182 deltas = [a for a in args if isinstance(a, KroneckerDelta)] 

1183 for delta1, delta2 in subsets(deltas, 2): 

1184 args1 = delta1.args 

1185 args2 = delta2.args 

1186 if args_cancel(args1, args2): 

1187 return S.Zero * m # In case of oo etc 

1188 return m 

1189 

1190 if not expr.has(KroneckerDelta): 

1191 return expr 

1192 

1193 if expr.has(Piecewise): 

1194 expr = expr.rewrite(KroneckerDelta) 

1195 

1196 newexpr = expr 

1197 expr = None 

1198 

1199 while newexpr != expr: 

1200 expr = newexpr 

1201 newexpr = expr.replace(lambda e: isinstance(e, Mul), cancel_kronecker_mul) 

1202 

1203 return expr 

1204 

1205 

1206def besselsimp(expr): 

1207 """ 

1208 Simplify bessel-type functions. 

1209 

1210 Explanation 

1211 =========== 

1212 

1213 This routine tries to simplify bessel-type functions. Currently it only 

1214 works on the Bessel J and I functions, however. It works by looking at all 

1215 such functions in turn, and eliminating factors of "I" and "-1" (actually 

1216 their polar equivalents) in front of the argument. Then, functions of 

1217 half-integer order are rewritten using strigonometric functions and 

1218 functions of integer order (> 1) are rewritten using functions 

1219 of low order. Finally, if the expression was changed, compute 

1220 factorization of the result with factor(). 

1221 

1222 >>> from sympy import besselj, besseli, besselsimp, polar_lift, I, S 

1223 >>> from sympy.abc import z, nu 

1224 >>> besselsimp(besselj(nu, z*polar_lift(-1))) 

1225 exp(I*pi*nu)*besselj(nu, z) 

1226 >>> besselsimp(besseli(nu, z*polar_lift(-I))) 

1227 exp(-I*pi*nu/2)*besselj(nu, z) 

1228 >>> besselsimp(besseli(S(-1)/2, z)) 

1229 sqrt(2)*cosh(z)/(sqrt(pi)*sqrt(z)) 

1230 >>> besselsimp(z*besseli(0, z) + z*(besseli(2, z))/2 + besseli(1, z)) 

1231 3*z*besseli(0, z)/2 

1232 """ 

1233 # TODO 

1234 # - better algorithm? 

1235 # - simplify (cos(pi*b)*besselj(b,z) - besselj(-b,z))/sin(pi*b) ... 

1236 # - use contiguity relations? 

1237 

1238 def replacer(fro, to, factors): 

1239 factors = set(factors) 

1240 

1241 def repl(nu, z): 

1242 if factors.intersection(Mul.make_args(z)): 

1243 return to(nu, z) 

1244 return fro(nu, z) 

1245 return repl 

1246 

1247 def torewrite(fro, to): 

1248 def tofunc(nu, z): 

1249 return fro(nu, z).rewrite(to) 

1250 return tofunc 

1251 

1252 def tominus(fro): 

1253 def tofunc(nu, z): 

1254 return exp(I*pi*nu)*fro(nu, exp_polar(-I*pi)*z) 

1255 return tofunc 

1256 

1257 orig_expr = expr 

1258 

1259 ifactors = [I, exp_polar(I*pi/2), exp_polar(-I*pi/2)] 

1260 expr = expr.replace( 

1261 besselj, replacer(besselj, 

1262 torewrite(besselj, besseli), ifactors)) 

1263 expr = expr.replace( 

1264 besseli, replacer(besseli, 

1265 torewrite(besseli, besselj), ifactors)) 

1266 

1267 minusfactors = [-1, exp_polar(I*pi)] 

1268 expr = expr.replace( 

1269 besselj, replacer(besselj, tominus(besselj), minusfactors)) 

1270 expr = expr.replace( 

1271 besseli, replacer(besseli, tominus(besseli), minusfactors)) 

1272 

1273 z0 = Dummy('z') 

1274 

1275 def expander(fro): 

1276 def repl(nu, z): 

1277 if (nu % 1) == S.Half: 

1278 return simplify(trigsimp(unpolarify( 

1279 fro(nu, z0).rewrite(besselj).rewrite(jn).expand( 

1280 func=True)).subs(z0, z))) 

1281 elif nu.is_Integer and nu > 1: 

1282 return fro(nu, z).expand(func=True) 

1283 return fro(nu, z) 

1284 return repl 

1285 

1286 expr = expr.replace(besselj, expander(besselj)) 

1287 expr = expr.replace(bessely, expander(bessely)) 

1288 expr = expr.replace(besseli, expander(besseli)) 

1289 expr = expr.replace(besselk, expander(besselk)) 

1290 

1291 def _bessel_simp_recursion(expr): 

1292 

1293 def _use_recursion(bessel, expr): 

1294 while True: 

1295 bessels = expr.find(lambda x: isinstance(x, bessel)) 

1296 try: 

1297 for ba in sorted(bessels, key=lambda x: re(x.args[0])): 

1298 a, x = ba.args 

1299 bap1 = bessel(a+1, x) 

1300 bap2 = bessel(a+2, x) 

1301 if expr.has(bap1) and expr.has(bap2): 

1302 expr = expr.subs(ba, 2*(a+1)/x*bap1 - bap2) 

1303 break 

1304 else: 

1305 return expr 

1306 except (ValueError, TypeError): 

1307 return expr 

1308 if expr.has(besselj): 

1309 expr = _use_recursion(besselj, expr) 

1310 if expr.has(bessely): 

1311 expr = _use_recursion(bessely, expr) 

1312 return expr 

1313 

1314 expr = _bessel_simp_recursion(expr) 

1315 if expr != orig_expr: 

1316 expr = expr.factor() 

1317 

1318 return expr 

1319 

1320 

1321def nthroot(expr, n, max_len=4, prec=15): 

1322 """ 

1323 Compute a real nth-root of a sum of surds. 

1324 

1325 Parameters 

1326 ========== 

1327 

1328 expr : sum of surds 

1329 n : integer 

1330 max_len : maximum number of surds passed as constants to ``nsimplify`` 

1331 

1332 Algorithm 

1333 ========= 

1334 

1335 First ``nsimplify`` is used to get a candidate root; if it is not a 

1336 root the minimal polynomial is computed; the answer is one of its 

1337 roots. 

1338 

1339 Examples 

1340 ======== 

1341 

1342 >>> from sympy.simplify.simplify import nthroot 

1343 >>> from sympy import sqrt 

1344 >>> nthroot(90 + 34*sqrt(7), 3) 

1345 sqrt(7) + 3 

1346 

1347 """ 

1348 expr = sympify(expr) 

1349 n = sympify(n) 

1350 p = expr**Rational(1, n) 

1351 if not n.is_integer: 

1352 return p 

1353 if not _is_sum_surds(expr): 

1354 return p 

1355 surds = [] 

1356 coeff_muls = [x.as_coeff_Mul() for x in expr.args] 

1357 for x, y in coeff_muls: 

1358 if not x.is_rational: 

1359 return p 

1360 if y is S.One: 

1361 continue 

1362 if not (y.is_Pow and y.exp == S.Half and y.base.is_integer): 

1363 return p 

1364 surds.append(y) 

1365 surds.sort() 

1366 surds = surds[:max_len] 

1367 if expr < 0 and n % 2 == 1: 

1368 p = (-expr)**Rational(1, n) 

1369 a = nsimplify(p, constants=surds) 

1370 res = a if _mexpand(a**n) == _mexpand(-expr) else p 

1371 return -res 

1372 a = nsimplify(p, constants=surds) 

1373 if _mexpand(a) is not _mexpand(p) and _mexpand(a**n) == _mexpand(expr): 

1374 return _mexpand(a) 

1375 expr = _nthroot_solve(expr, n, prec) 

1376 if expr is None: 

1377 return p 

1378 return expr 

1379 

1380 

1381def nsimplify(expr, constants=(), tolerance=None, full=False, rational=None, 

1382 rational_conversion='base10'): 

1383 """ 

1384 Find a simple representation for a number or, if there are free symbols or 

1385 if ``rational=True``, then replace Floats with their Rational equivalents. If 

1386 no change is made and rational is not False then Floats will at least be 

1387 converted to Rationals. 

1388 

1389 Explanation 

1390 =========== 

1391 

1392 For numerical expressions, a simple formula that numerically matches the 

1393 given numerical expression is sought (and the input should be possible 

1394 to evalf to a precision of at least 30 digits). 

1395 

1396 Optionally, a list of (rationally independent) constants to 

1397 include in the formula may be given. 

1398 

1399 A lower tolerance may be set to find less exact matches. If no tolerance 

1400 is given then the least precise value will set the tolerance (e.g. Floats 

1401 default to 15 digits of precision, so would be tolerance=10**-15). 

1402 

1403 With ``full=True``, a more extensive search is performed 

1404 (this is useful to find simpler numbers when the tolerance 

1405 is set low). 

1406 

1407 When converting to rational, if rational_conversion='base10' (the default), then 

1408 convert floats to rationals using their base-10 (string) representation. 

1409 When rational_conversion='exact' it uses the exact, base-2 representation. 

1410 

1411 Examples 

1412 ======== 

1413 

1414 >>> from sympy import nsimplify, sqrt, GoldenRatio, exp, I, pi 

1415 >>> nsimplify(4/(1+sqrt(5)), [GoldenRatio]) 

1416 -2 + 2*GoldenRatio 

1417 >>> nsimplify((1/(exp(3*pi*I/5)+1))) 

1418 1/2 - I*sqrt(sqrt(5)/10 + 1/4) 

1419 >>> nsimplify(I**I, [pi]) 

1420 exp(-pi/2) 

1421 >>> nsimplify(pi, tolerance=0.01) 

1422 22/7 

1423 

1424 >>> nsimplify(0.333333333333333, rational=True, rational_conversion='exact') 

1425 6004799503160655/18014398509481984 

1426 >>> nsimplify(0.333333333333333, rational=True) 

1427 1/3 

1428 

1429 See Also 

1430 ======== 

1431 

1432 sympy.core.function.nfloat 

1433 

1434 """ 

1435 try: 

1436 return sympify(as_int(expr)) 

1437 except (TypeError, ValueError): 

1438 pass 

1439 expr = sympify(expr).xreplace({ 

1440 Float('inf'): S.Infinity, 

1441 Float('-inf'): S.NegativeInfinity, 

1442 }) 

1443 if expr is S.Infinity or expr is S.NegativeInfinity: 

1444 return expr 

1445 if rational or expr.free_symbols: 

1446 return _real_to_rational(expr, tolerance, rational_conversion) 

1447 

1448 # SymPy's default tolerance for Rationals is 15; other numbers may have 

1449 # lower tolerances set, so use them to pick the largest tolerance if None 

1450 # was given 

1451 if tolerance is None: 

1452 tolerance = 10**-min([15] + 

1453 [mpmath.libmp.libmpf.prec_to_dps(n._prec) 

1454 for n in expr.atoms(Float)]) 

1455 # XXX should prec be set independent of tolerance or should it be computed 

1456 # from tolerance? 

1457 prec = 30 

1458 bprec = int(prec*3.33) 

1459 

1460 constants_dict = {} 

1461 for constant in constants: 

1462 constant = sympify(constant) 

1463 v = constant.evalf(prec) 

1464 if not v.is_Float: 

1465 raise ValueError("constants must be real-valued") 

1466 constants_dict[str(constant)] = v._to_mpmath(bprec) 

1467 

1468 exprval = expr.evalf(prec, chop=True) 

1469 re, im = exprval.as_real_imag() 

1470 

1471 # safety check to make sure that this evaluated to a number 

1472 if not (re.is_Number and im.is_Number): 

1473 return expr 

1474 

1475 def nsimplify_real(x): 

1476 orig = mpmath.mp.dps 

1477 xv = x._to_mpmath(bprec) 

1478 try: 

1479 # We'll be happy with low precision if a simple fraction 

1480 if not (tolerance or full): 

1481 mpmath.mp.dps = 15 

1482 rat = mpmath.pslq([xv, 1]) 

1483 if rat is not None: 

1484 return Rational(-int(rat[1]), int(rat[0])) 

1485 mpmath.mp.dps = prec 

1486 newexpr = mpmath.identify(xv, constants=constants_dict, 

1487 tol=tolerance, full=full) 

1488 if not newexpr: 

1489 raise ValueError 

1490 if full: 

1491 newexpr = newexpr[0] 

1492 expr = sympify(newexpr) 

1493 if x and not expr: # don't let x become 0 

1494 raise ValueError 

1495 if expr.is_finite is False and xv not in [mpmath.inf, mpmath.ninf]: 

1496 raise ValueError 

1497 return expr 

1498 finally: 

1499 # even though there are returns above, this is executed 

1500 # before leaving 

1501 mpmath.mp.dps = orig 

1502 try: 

1503 if re: 

1504 re = nsimplify_real(re) 

1505 if im: 

1506 im = nsimplify_real(im) 

1507 except ValueError: 

1508 if rational is None: 

1509 return _real_to_rational(expr, rational_conversion=rational_conversion) 

1510 return expr 

1511 

1512 rv = re + im*S.ImaginaryUnit 

1513 # if there was a change or rational is explicitly not wanted 

1514 # return the value, else return the Rational representation 

1515 if rv != expr or rational is False: 

1516 return rv 

1517 return _real_to_rational(expr, rational_conversion=rational_conversion) 

1518 

1519 

1520def _real_to_rational(expr, tolerance=None, rational_conversion='base10'): 

1521 """ 

1522 Replace all reals in expr with rationals. 

1523 

1524 Examples 

1525 ======== 

1526 

1527 >>> from sympy.simplify.simplify import _real_to_rational 

1528 >>> from sympy.abc import x 

1529 

1530 >>> _real_to_rational(.76 + .1*x**.5) 

1531 sqrt(x)/10 + 19/25 

1532 

1533 If rational_conversion='base10', this uses the base-10 string. If 

1534 rational_conversion='exact', the exact, base-2 representation is used. 

1535 

1536 >>> _real_to_rational(0.333333333333333, rational_conversion='exact') 

1537 6004799503160655/18014398509481984 

1538 >>> _real_to_rational(0.333333333333333) 

1539 1/3 

1540 

1541 """ 

1542 expr = _sympify(expr) 

1543 inf = Float('inf') 

1544 p = expr 

1545 reps = {} 

1546 reduce_num = None 

1547 if tolerance is not None and tolerance < 1: 

1548 reduce_num = ceiling(1/tolerance) 

1549 for fl in p.atoms(Float): 

1550 key = fl 

1551 if reduce_num is not None: 

1552 r = Rational(fl).limit_denominator(reduce_num) 

1553 elif (tolerance is not None and tolerance >= 1 and 

1554 fl.is_Integer is False): 

1555 r = Rational(tolerance*round(fl/tolerance) 

1556 ).limit_denominator(int(tolerance)) 

1557 else: 

1558 if rational_conversion == 'exact': 

1559 r = Rational(fl) 

1560 reps[key] = r 

1561 continue 

1562 elif rational_conversion != 'base10': 

1563 raise ValueError("rational_conversion must be 'base10' or 'exact'") 

1564 

1565 r = nsimplify(fl, rational=False) 

1566 # e.g. log(3).n() -> log(3) instead of a Rational 

1567 if fl and not r: 

1568 r = Rational(fl) 

1569 elif not r.is_Rational: 

1570 if fl in (inf, -inf): 

1571 r = S.ComplexInfinity 

1572 elif fl < 0: 

1573 fl = -fl 

1574 d = Pow(10, int(mpmath.log(fl)/mpmath.log(10))) 

1575 r = -Rational(str(fl/d))*d 

1576 elif fl > 0: 

1577 d = Pow(10, int(mpmath.log(fl)/mpmath.log(10))) 

1578 r = Rational(str(fl/d))*d 

1579 else: 

1580 r = S.Zero 

1581 reps[key] = r 

1582 return p.subs(reps, simultaneous=True) 

1583 

1584 

1585def clear_coefficients(expr, rhs=S.Zero): 

1586 """Return `p, r` where `p` is the expression obtained when Rational 

1587 additive and multiplicative coefficients of `expr` have been stripped 

1588 away in a naive fashion (i.e. without simplification). The operations 

1589 needed to remove the coefficients will be applied to `rhs` and returned 

1590 as `r`. 

1591 

1592 Examples 

1593 ======== 

1594 

1595 >>> from sympy.simplify.simplify import clear_coefficients 

1596 >>> from sympy.abc import x, y 

1597 >>> from sympy import Dummy 

1598 >>> expr = 4*y*(6*x + 3) 

1599 >>> clear_coefficients(expr - 2) 

1600 (y*(2*x + 1), 1/6) 

1601 

1602 When solving 2 or more expressions like `expr = a`, 

1603 `expr = b`, etc..., it is advantageous to provide a Dummy symbol 

1604 for `rhs` and simply replace it with `a`, `b`, etc... in `r`. 

1605 

1606 >>> rhs = Dummy('rhs') 

1607 >>> clear_coefficients(expr, rhs) 

1608 (y*(2*x + 1), _rhs/12) 

1609 >>> _[1].subs(rhs, 2) 

1610 1/6 

1611 """ 

1612 was = None 

1613 free = expr.free_symbols 

1614 if expr.is_Rational: 

1615 return (S.Zero, rhs - expr) 

1616 while expr and was != expr: 

1617 was = expr 

1618 m, expr = ( 

1619 expr.as_content_primitive() 

1620 if free else 

1621 factor_terms(expr).as_coeff_Mul(rational=True)) 

1622 rhs /= m 

1623 c, expr = expr.as_coeff_Add(rational=True) 

1624 rhs -= c 

1625 expr = signsimp(expr, evaluate = False) 

1626 if expr.could_extract_minus_sign(): 

1627 expr = -expr 

1628 rhs = -rhs 

1629 return expr, rhs 

1630 

1631def nc_simplify(expr, deep=True): 

1632 ''' 

1633 Simplify a non-commutative expression composed of multiplication 

1634 and raising to a power by grouping repeated subterms into one power. 

1635 Priority is given to simplifications that give the fewest number 

1636 of arguments in the end (for example, in a*b*a*b*c*a*b*c simplifying 

1637 to (a*b)**2*c*a*b*c gives 5 arguments while a*b*(a*b*c)**2 has 3). 

1638 If ``expr`` is a sum of such terms, the sum of the simplified terms 

1639 is returned. 

1640 

1641 Keyword argument ``deep`` controls whether or not subexpressions 

1642 nested deeper inside the main expression are simplified. See examples 

1643 below. Setting `deep` to `False` can save time on nested expressions 

1644 that do not need simplifying on all levels. 

1645 

1646 Examples 

1647 ======== 

1648 

1649 >>> from sympy import symbols 

1650 >>> from sympy.simplify.simplify import nc_simplify 

1651 >>> a, b, c = symbols("a b c", commutative=False) 

1652 >>> nc_simplify(a*b*a*b*c*a*b*c) 

1653 a*b*(a*b*c)**2 

1654 >>> expr = a**2*b*a**4*b*a**4 

1655 >>> nc_simplify(expr) 

1656 a**2*(b*a**4)**2 

1657 >>> nc_simplify(a*b*a*b*c**2*(a*b)**2*c**2) 

1658 ((a*b)**2*c**2)**2 

1659 >>> nc_simplify(a*b*a*b + 2*a*c*a**2*c*a**2*c*a) 

1660 (a*b)**2 + 2*(a*c*a)**3 

1661 >>> nc_simplify(b**-1*a**-1*(a*b)**2) 

1662 a*b 

1663 >>> nc_simplify(a**-1*b**-1*c*a) 

1664 (b*a)**(-1)*c*a 

1665 >>> expr = (a*b*a*b)**2*a*c*a*c 

1666 >>> nc_simplify(expr) 

1667 (a*b)**4*(a*c)**2 

1668 >>> nc_simplify(expr, deep=False) 

1669 (a*b*a*b)**2*(a*c)**2 

1670 

1671 ''' 

1672 if isinstance(expr, MatrixExpr): 

1673 expr = expr.doit(inv_expand=False) 

1674 _Add, _Mul, _Pow, _Symbol = MatAdd, MatMul, MatPow, MatrixSymbol 

1675 else: 

1676 _Add, _Mul, _Pow, _Symbol = Add, Mul, Pow, Symbol 

1677 

1678 # =========== Auxiliary functions ======================== 

1679 def _overlaps(args): 

1680 # Calculate a list of lists m such that m[i][j] contains the lengths 

1681 # of all possible overlaps between args[:i+1] and args[i+1+j:]. 

1682 # An overlap is a suffix of the prefix that matches a prefix 

1683 # of the suffix. 

1684 # For example, let expr=c*a*b*a*b*a*b*a*b. Then m[3][0] contains 

1685 # the lengths of overlaps of c*a*b*a*b with a*b*a*b. The overlaps 

1686 # are a*b*a*b, a*b and the empty word so that m[3][0]=[4,2,0]. 

1687 # All overlaps rather than only the longest one are recorded 

1688 # because this information helps calculate other overlap lengths. 

1689 m = [[([1, 0] if a == args[0] else [0]) for a in args[1:]]] 

1690 for i in range(1, len(args)): 

1691 overlaps = [] 

1692 j = 0 

1693 for j in range(len(args) - i - 1): 

1694 overlap = [] 

1695 for v in m[i-1][j+1]: 

1696 if j + i + 1 + v < len(args) and args[i] == args[j+i+1+v]: 

1697 overlap.append(v + 1) 

1698 overlap += [0] 

1699 overlaps.append(overlap) 

1700 m.append(overlaps) 

1701 return m 

1702 

1703 def _reduce_inverses(_args): 

1704 # replace consecutive negative powers by an inverse 

1705 # of a product of positive powers, e.g. a**-1*b**-1*c 

1706 # will simplify to (a*b)**-1*c; 

1707 # return that new args list and the number of negative 

1708 # powers in it (inv_tot) 

1709 inv_tot = 0 # total number of inverses 

1710 inverses = [] 

1711 args = [] 

1712 for arg in _args: 

1713 if isinstance(arg, _Pow) and arg.args[1].is_extended_negative: 

1714 inverses = [arg**-1] + inverses 

1715 inv_tot += 1 

1716 else: 

1717 if len(inverses) == 1: 

1718 args.append(inverses[0]**-1) 

1719 elif len(inverses) > 1: 

1720 args.append(_Pow(_Mul(*inverses), -1)) 

1721 inv_tot -= len(inverses) - 1 

1722 inverses = [] 

1723 args.append(arg) 

1724 if inverses: 

1725 args.append(_Pow(_Mul(*inverses), -1)) 

1726 inv_tot -= len(inverses) - 1 

1727 return inv_tot, tuple(args) 

1728 

1729 def get_score(s): 

1730 # compute the number of arguments of s 

1731 # (including in nested expressions) overall 

1732 # but ignore exponents 

1733 if isinstance(s, _Pow): 

1734 return get_score(s.args[0]) 

1735 elif isinstance(s, (_Add, _Mul)): 

1736 return sum([get_score(a) for a in s.args]) 

1737 return 1 

1738 

1739 def compare(s, alt_s): 

1740 # compare two possible simplifications and return a 

1741 # "better" one 

1742 if s != alt_s and get_score(alt_s) < get_score(s): 

1743 return alt_s 

1744 return s 

1745 # ======================================================== 

1746 

1747 if not isinstance(expr, (_Add, _Mul, _Pow)) or expr.is_commutative: 

1748 return expr 

1749 args = expr.args[:] 

1750 if isinstance(expr, _Pow): 

1751 if deep: 

1752 return _Pow(nc_simplify(args[0]), args[1]).doit() 

1753 else: 

1754 return expr 

1755 elif isinstance(expr, _Add): 

1756 return _Add(*[nc_simplify(a, deep=deep) for a in args]).doit() 

1757 else: 

1758 # get the non-commutative part 

1759 c_args, args = expr.args_cnc() 

1760 com_coeff = Mul(*c_args) 

1761 if com_coeff != 1: 

1762 return com_coeff*nc_simplify(expr/com_coeff, deep=deep) 

1763 

1764 inv_tot, args = _reduce_inverses(args) 

1765 # if most arguments are negative, work with the inverse 

1766 # of the expression, e.g. a**-1*b*a**-1*c**-1 will become 

1767 # (c*a*b**-1*a)**-1 at the end so can work with c*a*b**-1*a 

1768 invert = False 

1769 if inv_tot > len(args)/2: 

1770 invert = True 

1771 args = [a**-1 for a in args[::-1]] 

1772 

1773 if deep: 

1774 args = tuple(nc_simplify(a) for a in args) 

1775 

1776 m = _overlaps(args) 

1777 

1778 # simps will be {subterm: end} where `end` is the ending 

1779 # index of a sequence of repetitions of subterm; 

1780 # this is for not wasting time with subterms that are part 

1781 # of longer, already considered sequences 

1782 simps = {} 

1783 

1784 post = 1 

1785 pre = 1 

1786 

1787 # the simplification coefficient is the number of 

1788 # arguments by which contracting a given sequence 

1789 # would reduce the word; e.g. in a*b*a*b*c*a*b*c, 

1790 # contracting a*b*a*b to (a*b)**2 removes 3 arguments 

1791 # while a*b*c*a*b*c to (a*b*c)**2 removes 6. It's 

1792 # better to contract the latter so simplification 

1793 # with a maximum simplification coefficient will be chosen 

1794 max_simp_coeff = 0 

1795 simp = None # information about future simplification 

1796 

1797 for i in range(1, len(args)): 

1798 simp_coeff = 0 

1799 l = 0 # length of a subterm 

1800 p = 0 # the power of a subterm 

1801 if i < len(args) - 1: 

1802 rep = m[i][0] 

1803 start = i # starting index of the repeated sequence 

1804 end = i+1 # ending index of the repeated sequence 

1805 if i == len(args)-1 or rep == [0]: 

1806 # no subterm is repeated at this stage, at least as 

1807 # far as the arguments are concerned - there may be 

1808 # a repetition if powers are taken into account 

1809 if (isinstance(args[i], _Pow) and 

1810 not isinstance(args[i].args[0], _Symbol)): 

1811 subterm = args[i].args[0].args 

1812 l = len(subterm) 

1813 if args[i-l:i] == subterm: 

1814 # e.g. a*b in a*b*(a*b)**2 is not repeated 

1815 # in args (= [a, b, (a*b)**2]) but it 

1816 # can be matched here 

1817 p += 1 

1818 start -= l 

1819 if args[i+1:i+1+l] == subterm: 

1820 # e.g. a*b in (a*b)**2*a*b 

1821 p += 1 

1822 end += l 

1823 if p: 

1824 p += args[i].args[1] 

1825 else: 

1826 continue 

1827 else: 

1828 l = rep[0] # length of the longest repeated subterm at this point 

1829 start -= l - 1 

1830 subterm = args[start:end] 

1831 p = 2 

1832 end += l 

1833 

1834 if subterm in simps and simps[subterm] >= start: 

1835 # the subterm is part of a sequence that 

1836 # has already been considered 

1837 continue 

1838 

1839 # count how many times it's repeated 

1840 while end < len(args): 

1841 if l in m[end-1][0]: 

1842 p += 1 

1843 end += l 

1844 elif isinstance(args[end], _Pow) and args[end].args[0].args == subterm: 

1845 # for cases like a*b*a*b*(a*b)**2*a*b 

1846 p += args[end].args[1] 

1847 end += 1 

1848 else: 

1849 break 

1850 

1851 # see if another match can be made, e.g. 

1852 # for b*a**2 in b*a**2*b*a**3 or a*b in 

1853 # a**2*b*a*b 

1854 

1855 pre_exp = 0 

1856 pre_arg = 1 

1857 if start - l >= 0 and args[start-l+1:start] == subterm[1:]: 

1858 if isinstance(subterm[0], _Pow): 

1859 pre_arg = subterm[0].args[0] 

1860 exp = subterm[0].args[1] 

1861 else: 

1862 pre_arg = subterm[0] 

1863 exp = 1 

1864 if isinstance(args[start-l], _Pow) and args[start-l].args[0] == pre_arg: 

1865 pre_exp = args[start-l].args[1] - exp 

1866 start -= l 

1867 p += 1 

1868 elif args[start-l] == pre_arg: 

1869 pre_exp = 1 - exp 

1870 start -= l 

1871 p += 1 

1872 

1873 post_exp = 0 

1874 post_arg = 1 

1875 if end + l - 1 < len(args) and args[end:end+l-1] == subterm[:-1]: 

1876 if isinstance(subterm[-1], _Pow): 

1877 post_arg = subterm[-1].args[0] 

1878 exp = subterm[-1].args[1] 

1879 else: 

1880 post_arg = subterm[-1] 

1881 exp = 1 

1882 if isinstance(args[end+l-1], _Pow) and args[end+l-1].args[0] == post_arg: 

1883 post_exp = args[end+l-1].args[1] - exp 

1884 end += l 

1885 p += 1 

1886 elif args[end+l-1] == post_arg: 

1887 post_exp = 1 - exp 

1888 end += l 

1889 p += 1 

1890 

1891 # Consider a*b*a**2*b*a**2*b*a: 

1892 # b*a**2 is explicitly repeated, but note 

1893 # that in this case a*b*a is also repeated 

1894 # so there are two possible simplifications: 

1895 # a*(b*a**2)**3*a**-1 or (a*b*a)**3 

1896 # The latter is obviously simpler. 

1897 # But in a*b*a**2*b**2*a**2 the simplifications are 

1898 # a*(b*a**2)**2 and (a*b*a)**3*a in which case 

1899 # it's better to stick with the shorter subterm 

1900 if post_exp and exp % 2 == 0 and start > 0: 

1901 exp = exp/2 

1902 _pre_exp = 1 

1903 _post_exp = 1 

1904 if isinstance(args[start-1], _Pow) and args[start-1].args[0] == post_arg: 

1905 _post_exp = post_exp + exp 

1906 _pre_exp = args[start-1].args[1] - exp 

1907 elif args[start-1] == post_arg: 

1908 _post_exp = post_exp + exp 

1909 _pre_exp = 1 - exp 

1910 if _pre_exp == 0 or _post_exp == 0: 

1911 if not pre_exp: 

1912 start -= 1 

1913 post_exp = _post_exp 

1914 pre_exp = _pre_exp 

1915 pre_arg = post_arg 

1916 subterm = (post_arg**exp,) + subterm[:-1] + (post_arg**exp,) 

1917 

1918 simp_coeff += end-start 

1919 

1920 if post_exp: 

1921 simp_coeff -= 1 

1922 if pre_exp: 

1923 simp_coeff -= 1 

1924 

1925 simps[subterm] = end 

1926 

1927 if simp_coeff > max_simp_coeff: 

1928 max_simp_coeff = simp_coeff 

1929 simp = (start, _Mul(*subterm), p, end, l) 

1930 pre = pre_arg**pre_exp 

1931 post = post_arg**post_exp 

1932 

1933 if simp: 

1934 subterm = _Pow(nc_simplify(simp[1], deep=deep), simp[2]) 

1935 pre = nc_simplify(_Mul(*args[:simp[0]])*pre, deep=deep) 

1936 post = post*nc_simplify(_Mul(*args[simp[3]:]), deep=deep) 

1937 simp = pre*subterm*post 

1938 if pre != 1 or post != 1: 

1939 # new simplifications may be possible but no need 

1940 # to recurse over arguments 

1941 simp = nc_simplify(simp, deep=False) 

1942 else: 

1943 simp = _Mul(*args) 

1944 

1945 if invert: 

1946 simp = _Pow(simp, -1) 

1947 

1948 # see if factor_nc(expr) is simplified better 

1949 if not isinstance(expr, MatrixExpr): 

1950 f_expr = factor_nc(expr) 

1951 if f_expr != expr: 

1952 alt_simp = nc_simplify(f_expr, deep=deep) 

1953 simp = compare(simp, alt_simp) 

1954 else: 

1955 simp = simp.doit(inv_expand=False) 

1956 return simp 

1957 

1958 

1959def dotprodsimp(expr, withsimp=False): 

1960 """Simplification for a sum of products targeted at the kind of blowup that 

1961 occurs during summation of products. Intended to reduce expression blowup 

1962 during matrix multiplication or other similar operations. Only works with 

1963 algebraic expressions and does not recurse into non. 

1964 

1965 Parameters 

1966 ========== 

1967 

1968 withsimp : bool, optional 

1969 Specifies whether a flag should be returned along with the expression 

1970 to indicate roughly whether simplification was successful. It is used 

1971 in ``MatrixArithmetic._eval_pow_by_recursion`` to avoid attempting to 

1972 simplify an expression repetitively which does not simplify. 

1973 """ 

1974 

1975 def count_ops_alg(expr): 

1976 """Optimized count algebraic operations with no recursion into 

1977 non-algebraic args that ``core.function.count_ops`` does. Also returns 

1978 whether rational functions may be present according to negative 

1979 exponents of powers or non-number fractions. 

1980 

1981 Returns 

1982 ======= 

1983 

1984 ops, ratfunc : int, bool 

1985 ``ops`` is the number of algebraic operations starting at the top 

1986 level expression (not recursing into non-alg children). ``ratfunc`` 

1987 specifies whether the expression MAY contain rational functions 

1988 which ``cancel`` MIGHT optimize. 

1989 """ 

1990 

1991 ops = 0 

1992 args = [expr] 

1993 ratfunc = False 

1994 

1995 while args: 

1996 a = args.pop() 

1997 

1998 if not isinstance(a, Basic): 

1999 continue 

2000 

2001 if a.is_Rational: 

2002 if a is not S.One: # -1/3 = NEG + DIV 

2003 ops += bool (a.p < 0) + bool (a.q != 1) 

2004 

2005 elif a.is_Mul: 

2006 if a.could_extract_minus_sign(): 

2007 ops += 1 

2008 if a.args[0] is S.NegativeOne: 

2009 a = a.as_two_terms()[1] 

2010 else: 

2011 a = -a 

2012 

2013 n, d = fraction(a) 

2014 

2015 if n.is_Integer: 

2016 ops += 1 + bool (n < 0) 

2017 args.append(d) # won't be -Mul but could be Add 

2018 

2019 elif d is not S.One: 

2020 if not d.is_Integer: 

2021 args.append(d) 

2022 ratfunc=True 

2023 

2024 ops += 1 

2025 args.append(n) # could be -Mul 

2026 

2027 else: 

2028 ops += len(a.args) - 1 

2029 args.extend(a.args) 

2030 

2031 elif a.is_Add: 

2032 laargs = len(a.args) 

2033 negs = 0 

2034 

2035 for ai in a.args: 

2036 if ai.could_extract_minus_sign(): 

2037 negs += 1 

2038 ai = -ai 

2039 args.append(ai) 

2040 

2041 ops += laargs - (negs != laargs) # -x - y = NEG + SUB 

2042 

2043 elif a.is_Pow: 

2044 ops += 1 

2045 args.append(a.base) 

2046 

2047 if not ratfunc: 

2048 ratfunc = a.exp.is_negative is not False 

2049 

2050 return ops, ratfunc 

2051 

2052 def nonalg_subs_dummies(expr, dummies): 

2053 """Substitute dummy variables for non-algebraic expressions to avoid 

2054 evaluation of non-algebraic terms that ``polys.polytools.cancel`` does. 

2055 """ 

2056 

2057 if not expr.args: 

2058 return expr 

2059 

2060 if expr.is_Add or expr.is_Mul or expr.is_Pow: 

2061 args = None 

2062 

2063 for i, a in enumerate(expr.args): 

2064 c = nonalg_subs_dummies(a, dummies) 

2065 

2066 if c is a: 

2067 continue 

2068 

2069 if args is None: 

2070 args = list(expr.args) 

2071 

2072 args[i] = c 

2073 

2074 if args is None: 

2075 return expr 

2076 

2077 return expr.func(*args) 

2078 

2079 return dummies.setdefault(expr, Dummy()) 

2080 

2081 simplified = False # doesn't really mean simplified, rather "can simplify again" 

2082 

2083 if isinstance(expr, Basic) and (expr.is_Add or expr.is_Mul or expr.is_Pow): 

2084 expr2 = expr.expand(deep=True, modulus=None, power_base=False, 

2085 power_exp=False, mul=True, log=False, multinomial=True, basic=False) 

2086 

2087 if expr2 != expr: 

2088 expr = expr2 

2089 simplified = True 

2090 

2091 exprops, ratfunc = count_ops_alg(expr) 

2092 

2093 if exprops >= 6: # empirically tested cutoff for expensive simplification 

2094 if ratfunc: 

2095 dummies = {} 

2096 expr2 = nonalg_subs_dummies(expr, dummies) 

2097 

2098 if expr2 is expr or count_ops_alg(expr2)[0] >= 6: # check again after substitution 

2099 expr3 = cancel(expr2) 

2100 

2101 if expr3 != expr2: 

2102 expr = expr3.subs([(d, e) for e, d in dummies.items()]) 

2103 simplified = True 

2104 

2105 # very special case: x/(x-1) - 1/(x-1) -> 1 

2106 elif (exprops == 5 and expr.is_Add and expr.args [0].is_Mul and 

2107 expr.args [1].is_Mul and expr.args [0].args [-1].is_Pow and 

2108 expr.args [1].args [-1].is_Pow and 

2109 expr.args [0].args [-1].exp is S.NegativeOne and 

2110 expr.args [1].args [-1].exp is S.NegativeOne): 

2111 

2112 expr2 = together (expr) 

2113 expr2ops = count_ops_alg(expr2)[0] 

2114 

2115 if expr2ops < exprops: 

2116 expr = expr2 

2117 simplified = True 

2118 

2119 else: 

2120 simplified = True 

2121 

2122 return (expr, simplified) if withsimp else expr 

2123 

2124 

2125bottom_up = deprecated( 

2126 """ 

2127 Using bottom_up from the sympy.simplify.simplify submodule is 

2128 deprecated. 

2129 

2130 Instead, use bottom_up from the top-level sympy namespace, like 

2131 

2132 sympy.bottom_up 

2133 """, 

2134 deprecated_since_version="1.10", 

2135 active_deprecations_target="deprecated-traversal-functions-moved", 

2136)(_bottom_up) 

2137 

2138 

2139# XXX: This function really should either be private API or exported in the 

2140# top-level sympy/__init__.py 

2141walk = deprecated( 

2142 """ 

2143 Using walk from the sympy.simplify.simplify submodule is 

2144 deprecated. 

2145 

2146 Instead, use walk from sympy.core.traversal.walk 

2147 """, 

2148 deprecated_since_version="1.10", 

2149 active_deprecations_target="deprecated-traversal-functions-moved", 

2150)(_walk)