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

512 statements  

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

1from collections import defaultdict 

2 

3from sympy.core import sympify, S, Mul, Derivative, Pow 

4from sympy.core.add import _unevaluated_Add, Add 

5from sympy.core.assumptions import assumptions 

6from sympy.core.exprtools import Factors, gcd_terms 

7from sympy.core.function import _mexpand, expand_mul, expand_power_base 

8from sympy.core.mul import _keep_coeff, _unevaluated_Mul, _mulsort 

9from sympy.core.numbers import Rational, zoo, nan 

10from sympy.core.parameters import global_parameters 

11from sympy.core.sorting import ordered, default_sort_key 

12from sympy.core.symbol import Dummy, Wild, symbols 

13from sympy.functions import exp, sqrt, log 

14from sympy.functions.elementary.complexes import Abs 

15from sympy.polys import gcd 

16from sympy.simplify.sqrtdenest import sqrtdenest 

17from sympy.utilities.iterables import iterable, sift 

18 

19 

20 

21 

22def collect(expr, syms, func=None, evaluate=None, exact=False, distribute_order_term=True): 

23 """ 

24 Collect additive terms of an expression. 

25 

26 Explanation 

27 =========== 

28 

29 This function collects additive terms of an expression with respect 

30 to a list of expression up to powers with rational exponents. By the 

31 term symbol here are meant arbitrary expressions, which can contain 

32 powers, products, sums etc. In other words symbol is a pattern which 

33 will be searched for in the expression's terms. 

34 

35 The input expression is not expanded by :func:`collect`, so user is 

36 expected to provide an expression in an appropriate form. This makes 

37 :func:`collect` more predictable as there is no magic happening behind the 

38 scenes. However, it is important to note, that powers of products are 

39 converted to products of powers using the :func:`~.expand_power_base` 

40 function. 

41 

42 There are two possible types of output. First, if ``evaluate`` flag is 

43 set, this function will return an expression with collected terms or 

44 else it will return a dictionary with expressions up to rational powers 

45 as keys and collected coefficients as values. 

46 

47 Examples 

48 ======== 

49 

50 >>> from sympy import S, collect, expand, factor, Wild 

51 >>> from sympy.abc import a, b, c, x, y 

52 

53 This function can collect symbolic coefficients in polynomials or 

54 rational expressions. It will manage to find all integer or rational 

55 powers of collection variable:: 

56 

57 >>> collect(a*x**2 + b*x**2 + a*x - b*x + c, x) 

58 c + x**2*(a + b) + x*(a - b) 

59 

60 The same result can be achieved in dictionary form:: 

61 

62 >>> d = collect(a*x**2 + b*x**2 + a*x - b*x + c, x, evaluate=False) 

63 >>> d[x**2] 

64 a + b 

65 >>> d[x] 

66 a - b 

67 >>> d[S.One] 

68 c 

69 

70 You can also work with multivariate polynomials. However, remember that 

71 this function is greedy so it will care only about a single symbol at time, 

72 in specification order:: 

73 

74 >>> collect(x**2 + y*x**2 + x*y + y + a*y, [x, y]) 

75 x**2*(y + 1) + x*y + y*(a + 1) 

76 

77 Also more complicated expressions can be used as patterns:: 

78 

79 >>> from sympy import sin, log 

80 >>> collect(a*sin(2*x) + b*sin(2*x), sin(2*x)) 

81 (a + b)*sin(2*x) 

82 

83 >>> collect(a*x*log(x) + b*(x*log(x)), x*log(x)) 

84 x*(a + b)*log(x) 

85 

86 You can use wildcards in the pattern:: 

87 

88 >>> w = Wild('w1') 

89 >>> collect(a*x**y - b*x**y, w**y) 

90 x**y*(a - b) 

91 

92 It is also possible to work with symbolic powers, although it has more 

93 complicated behavior, because in this case power's base and symbolic part 

94 of the exponent are treated as a single symbol:: 

95 

96 >>> collect(a*x**c + b*x**c, x) 

97 a*x**c + b*x**c 

98 >>> collect(a*x**c + b*x**c, x**c) 

99 x**c*(a + b) 

100 

101 However if you incorporate rationals to the exponents, then you will get 

102 well known behavior:: 

103 

104 >>> collect(a*x**(2*c) + b*x**(2*c), x**c) 

105 x**(2*c)*(a + b) 

106 

107 Note also that all previously stated facts about :func:`collect` function 

108 apply to the exponential function, so you can get:: 

109 

110 >>> from sympy import exp 

111 >>> collect(a*exp(2*x) + b*exp(2*x), exp(x)) 

112 (a + b)*exp(2*x) 

113 

114 If you are interested only in collecting specific powers of some symbols 

115 then set ``exact`` flag to True:: 

116 

117 >>> collect(a*x**7 + b*x**7, x, exact=True) 

118 a*x**7 + b*x**7 

119 >>> collect(a*x**7 + b*x**7, x**7, exact=True) 

120 x**7*(a + b) 

121 

122 If you want to collect on any object containing symbols, set 

123 ``exact`` to None: 

124 

125 >>> collect(x*exp(x) + sin(x)*y + sin(x)*2 + 3*x, x, exact=None) 

126 x*exp(x) + 3*x + (y + 2)*sin(x) 

127 >>> collect(a*x*y + x*y + b*x + x, [x, y], exact=None) 

128 x*y*(a + 1) + x*(b + 1) 

129 

130 You can also apply this function to differential equations, where 

131 derivatives of arbitrary order can be collected. Note that if you 

132 collect with respect to a function or a derivative of a function, all 

133 derivatives of that function will also be collected. Use 

134 ``exact=True`` to prevent this from happening:: 

135 

136 >>> from sympy import Derivative as D, collect, Function 

137 >>> f = Function('f') (x) 

138 

139 >>> collect(a*D(f,x) + b*D(f,x), D(f,x)) 

140 (a + b)*Derivative(f(x), x) 

141 

142 >>> collect(a*D(D(f,x),x) + b*D(D(f,x),x), f) 

143 (a + b)*Derivative(f(x), (x, 2)) 

144 

145 >>> collect(a*D(D(f,x),x) + b*D(D(f,x),x), D(f,x), exact=True) 

146 a*Derivative(f(x), (x, 2)) + b*Derivative(f(x), (x, 2)) 

147 

148 >>> collect(a*D(f,x) + b*D(f,x) + a*f + b*f, f) 

149 (a + b)*f(x) + (a + b)*Derivative(f(x), x) 

150 

151 Or you can even match both derivative order and exponent at the same time:: 

152 

153 >>> collect(a*D(D(f,x),x)**2 + b*D(D(f,x),x)**2, D(f,x)) 

154 (a + b)*Derivative(f(x), (x, 2))**2 

155 

156 Finally, you can apply a function to each of the collected coefficients. 

157 For example you can factorize symbolic coefficients of polynomial:: 

158 

159 >>> f = expand((x + a + 1)**3) 

160 

161 >>> collect(f, x, factor) 

162 x**3 + 3*x**2*(a + 1) + 3*x*(a + 1)**2 + (a + 1)**3 

163 

164 .. note:: Arguments are expected to be in expanded form, so you might have 

165 to call :func:`~.expand` prior to calling this function. 

166 

167 See Also 

168 ======== 

169 

170 collect_const, collect_sqrt, rcollect 

171 """ 

172 expr = sympify(expr) 

173 syms = [sympify(i) for i in (syms if iterable(syms) else [syms])] 

174 

175 # replace syms[i] if it is not x, -x or has Wild symbols 

176 cond = lambda x: x.is_Symbol or (-x).is_Symbol or bool( 

177 x.atoms(Wild)) 

178 _, nonsyms = sift(syms, cond, binary=True) 

179 if nonsyms: 

180 reps = dict(zip(nonsyms, [Dummy(**assumptions(i)) for i in nonsyms])) 

181 syms = [reps.get(s, s) for s in syms] 

182 rv = collect(expr.subs(reps), syms, 

183 func=func, evaluate=evaluate, exact=exact, 

184 distribute_order_term=distribute_order_term) 

185 urep = {v: k for k, v in reps.items()} 

186 if not isinstance(rv, dict): 

187 return rv.xreplace(urep) 

188 else: 

189 return {urep.get(k, k).xreplace(urep): v.xreplace(urep) 

190 for k, v in rv.items()} 

191 

192 # see if other expressions should be considered 

193 if exact is None: 

194 _syms = set() 

195 for i in Add.make_args(expr): 

196 if not i.has_free(*syms) or i in syms: 

197 continue 

198 if not i.is_Mul and i not in syms: 

199 _syms.add(i) 

200 else: 

201 # identify compound generators 

202 g = i._new_rawargs(*i.as_coeff_mul(*syms)[1]) 

203 if g not in syms: 

204 _syms.add(g) 

205 simple = all(i.is_Pow and i.base in syms for i in _syms) 

206 syms = syms + list(ordered(_syms)) 

207 if not simple: 

208 return collect(expr, syms, 

209 func=func, evaluate=evaluate, exact=False, 

210 distribute_order_term=distribute_order_term) 

211 

212 if evaluate is None: 

213 evaluate = global_parameters.evaluate 

214 

215 def make_expression(terms): 

216 product = [] 

217 

218 for term, rat, sym, deriv in terms: 

219 if deriv is not None: 

220 var, order = deriv 

221 

222 while order > 0: 

223 term, order = Derivative(term, var), order - 1 

224 

225 if sym is None: 

226 if rat is S.One: 

227 product.append(term) 

228 else: 

229 product.append(Pow(term, rat)) 

230 else: 

231 product.append(Pow(term, rat*sym)) 

232 

233 return Mul(*product) 

234 

235 def parse_derivative(deriv): 

236 # scan derivatives tower in the input expression and return 

237 # underlying function and maximal differentiation order 

238 expr, sym, order = deriv.expr, deriv.variables[0], 1 

239 

240 for s in deriv.variables[1:]: 

241 if s == sym: 

242 order += 1 

243 else: 

244 raise NotImplementedError( 

245 'Improve MV Derivative support in collect') 

246 

247 while isinstance(expr, Derivative): 

248 s0 = expr.variables[0] 

249 

250 for s in expr.variables: 

251 if s != s0: 

252 raise NotImplementedError( 

253 'Improve MV Derivative support in collect') 

254 

255 if s0 == sym: 

256 expr, order = expr.expr, order + len(expr.variables) 

257 else: 

258 break 

259 

260 return expr, (sym, Rational(order)) 

261 

262 def parse_term(expr): 

263 """Parses expression expr and outputs tuple (sexpr, rat_expo, 

264 sym_expo, deriv) 

265 where: 

266 - sexpr is the base expression 

267 - rat_expo is the rational exponent that sexpr is raised to 

268 - sym_expo is the symbolic exponent that sexpr is raised to 

269 - deriv contains the derivatives of the expression 

270 

271 For example, the output of x would be (x, 1, None, None) 

272 the output of 2**x would be (2, 1, x, None). 

273 """ 

274 rat_expo, sym_expo = S.One, None 

275 sexpr, deriv = expr, None 

276 

277 if expr.is_Pow: 

278 if isinstance(expr.base, Derivative): 

279 sexpr, deriv = parse_derivative(expr.base) 

280 else: 

281 sexpr = expr.base 

282 

283 if expr.base == S.Exp1: 

284 arg = expr.exp 

285 if arg.is_Rational: 

286 sexpr, rat_expo = S.Exp1, arg 

287 elif arg.is_Mul: 

288 coeff, tail = arg.as_coeff_Mul(rational=True) 

289 sexpr, rat_expo = exp(tail), coeff 

290 

291 elif expr.exp.is_Number: 

292 rat_expo = expr.exp 

293 else: 

294 coeff, tail = expr.exp.as_coeff_Mul() 

295 

296 if coeff.is_Number: 

297 rat_expo, sym_expo = coeff, tail 

298 else: 

299 sym_expo = expr.exp 

300 elif isinstance(expr, exp): 

301 arg = expr.exp 

302 if arg.is_Rational: 

303 sexpr, rat_expo = S.Exp1, arg 

304 elif arg.is_Mul: 

305 coeff, tail = arg.as_coeff_Mul(rational=True) 

306 sexpr, rat_expo = exp(tail), coeff 

307 elif isinstance(expr, Derivative): 

308 sexpr, deriv = parse_derivative(expr) 

309 

310 return sexpr, rat_expo, sym_expo, deriv 

311 

312 def parse_expression(terms, pattern): 

313 """Parse terms searching for a pattern. 

314 Terms is a list of tuples as returned by parse_terms; 

315 Pattern is an expression treated as a product of factors. 

316 """ 

317 pattern = Mul.make_args(pattern) 

318 

319 if len(terms) < len(pattern): 

320 # pattern is longer than matched product 

321 # so no chance for positive parsing result 

322 return None 

323 else: 

324 pattern = [parse_term(elem) for elem in pattern] 

325 

326 terms = terms[:] # need a copy 

327 elems, common_expo, has_deriv = [], None, False 

328 

329 for elem, e_rat, e_sym, e_ord in pattern: 

330 

331 if elem.is_Number and e_rat == 1 and e_sym is None: 

332 # a constant is a match for everything 

333 continue 

334 

335 for j in range(len(terms)): 

336 if terms[j] is None: 

337 continue 

338 

339 term, t_rat, t_sym, t_ord = terms[j] 

340 

341 # keeping track of whether one of the terms had 

342 # a derivative or not as this will require rebuilding 

343 # the expression later 

344 if t_ord is not None: 

345 has_deriv = True 

346 

347 if (term.match(elem) is not None and 

348 (t_sym == e_sym or t_sym is not None and 

349 e_sym is not None and 

350 t_sym.match(e_sym) is not None)): 

351 if exact is False: 

352 # we don't have to be exact so find common exponent 

353 # for both expression's term and pattern's element 

354 expo = t_rat / e_rat 

355 

356 if common_expo is None: 

357 # first time 

358 common_expo = expo 

359 else: 

360 # common exponent was negotiated before so 

361 # there is no chance for a pattern match unless 

362 # common and current exponents are equal 

363 if common_expo != expo: 

364 common_expo = 1 

365 else: 

366 # we ought to be exact so all fields of 

367 # interest must match in every details 

368 if e_rat != t_rat or e_ord != t_ord: 

369 continue 

370 

371 # found common term so remove it from the expression 

372 # and try to match next element in the pattern 

373 elems.append(terms[j]) 

374 terms[j] = None 

375 

376 break 

377 

378 else: 

379 # pattern element not found 

380 return None 

381 

382 return [_f for _f in terms if _f], elems, common_expo, has_deriv 

383 

384 if evaluate: 

385 if expr.is_Add: 

386 o = expr.getO() or 0 

387 expr = expr.func(*[ 

388 collect(a, syms, func, True, exact, distribute_order_term) 

389 for a in expr.args if a != o]) + o 

390 elif expr.is_Mul: 

391 return expr.func(*[ 

392 collect(term, syms, func, True, exact, distribute_order_term) 

393 for term in expr.args]) 

394 elif expr.is_Pow: 

395 b = collect( 

396 expr.base, syms, func, True, exact, distribute_order_term) 

397 return Pow(b, expr.exp) 

398 

399 syms = [expand_power_base(i, deep=False) for i in syms] 

400 

401 order_term = None 

402 

403 if distribute_order_term: 

404 order_term = expr.getO() 

405 

406 if order_term is not None: 

407 if order_term.has(*syms): 

408 order_term = None 

409 else: 

410 expr = expr.removeO() 

411 

412 summa = [expand_power_base(i, deep=False) for i in Add.make_args(expr)] 

413 

414 collected, disliked = defaultdict(list), S.Zero 

415 for product in summa: 

416 c, nc = product.args_cnc(split_1=False) 

417 args = list(ordered(c)) + nc 

418 terms = [parse_term(i) for i in args] 

419 small_first = True 

420 

421 for symbol in syms: 

422 if isinstance(symbol, Derivative) and small_first: 

423 terms = list(reversed(terms)) 

424 small_first = not small_first 

425 result = parse_expression(terms, symbol) 

426 

427 if result is not None: 

428 if not symbol.is_commutative: 

429 raise AttributeError("Can not collect noncommutative symbol") 

430 

431 terms, elems, common_expo, has_deriv = result 

432 

433 # when there was derivative in current pattern we 

434 # will need to rebuild its expression from scratch 

435 if not has_deriv: 

436 margs = [] 

437 for elem in elems: 

438 if elem[2] is None: 

439 e = elem[1] 

440 else: 

441 e = elem[1]*elem[2] 

442 margs.append(Pow(elem[0], e)) 

443 index = Mul(*margs) 

444 else: 

445 index = make_expression(elems) 

446 terms = expand_power_base(make_expression(terms), deep=False) 

447 index = expand_power_base(index, deep=False) 

448 collected[index].append(terms) 

449 break 

450 else: 

451 # none of the patterns matched 

452 disliked += product 

453 # add terms now for each key 

454 collected = {k: Add(*v) for k, v in collected.items()} 

455 

456 if disliked is not S.Zero: 

457 collected[S.One] = disliked 

458 

459 if order_term is not None: 

460 for key, val in collected.items(): 

461 collected[key] = val + order_term 

462 

463 if func is not None: 

464 collected = { 

465 key: func(val) for key, val in collected.items()} 

466 

467 if evaluate: 

468 return Add(*[key*val for key, val in collected.items()]) 

469 else: 

470 return collected 

471 

472 

473def rcollect(expr, *vars): 

474 """ 

475 Recursively collect sums in an expression. 

476 

477 Examples 

478 ======== 

479 

480 >>> from sympy.simplify import rcollect 

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

482 

483 >>> expr = (x**2*y + x*y + x + y)/(x + y) 

484 

485 >>> rcollect(expr, y) 

486 (x + y*(x**2 + x + 1))/(x + y) 

487 

488 See Also 

489 ======== 

490 

491 collect, collect_const, collect_sqrt 

492 """ 

493 if expr.is_Atom or not expr.has(*vars): 

494 return expr 

495 else: 

496 expr = expr.__class__(*[rcollect(arg, *vars) for arg in expr.args]) 

497 

498 if expr.is_Add: 

499 return collect(expr, vars) 

500 else: 

501 return expr 

502 

503 

504def collect_sqrt(expr, evaluate=None): 

505 """Return expr with terms having common square roots collected together. 

506 If ``evaluate`` is False a count indicating the number of sqrt-containing 

507 terms will be returned and, if non-zero, the terms of the Add will be 

508 returned, else the expression itself will be returned as a single term. 

509 If ``evaluate`` is True, the expression with any collected terms will be 

510 returned. 

511 

512 Note: since I = sqrt(-1), it is collected, too. 

513 

514 Examples 

515 ======== 

516 

517 >>> from sympy import sqrt 

518 >>> from sympy.simplify.radsimp import collect_sqrt 

519 >>> from sympy.abc import a, b 

520 

521 >>> r2, r3, r5 = [sqrt(i) for i in [2, 3, 5]] 

522 >>> collect_sqrt(a*r2 + b*r2) 

523 sqrt(2)*(a + b) 

524 >>> collect_sqrt(a*r2 + b*r2 + a*r3 + b*r3) 

525 sqrt(2)*(a + b) + sqrt(3)*(a + b) 

526 >>> collect_sqrt(a*r2 + b*r2 + a*r3 + b*r5) 

527 sqrt(3)*a + sqrt(5)*b + sqrt(2)*(a + b) 

528 

529 If evaluate is False then the arguments will be sorted and 

530 returned as a list and a count of the number of sqrt-containing 

531 terms will be returned: 

532 

533 >>> collect_sqrt(a*r2 + b*r2 + a*r3 + b*r5, evaluate=False) 

534 ((sqrt(3)*a, sqrt(5)*b, sqrt(2)*(a + b)), 3) 

535 >>> collect_sqrt(a*sqrt(2) + b, evaluate=False) 

536 ((b, sqrt(2)*a), 1) 

537 >>> collect_sqrt(a + b, evaluate=False) 

538 ((a + b,), 0) 

539 

540 See Also 

541 ======== 

542 

543 collect, collect_const, rcollect 

544 """ 

545 if evaluate is None: 

546 evaluate = global_parameters.evaluate 

547 # this step will help to standardize any complex arguments 

548 # of sqrts 

549 coeff, expr = expr.as_content_primitive() 

550 vars = set() 

551 for a in Add.make_args(expr): 

552 for m in a.args_cnc()[0]: 

553 if m.is_number and ( 

554 m.is_Pow and m.exp.is_Rational and m.exp.q == 2 or 

555 m is S.ImaginaryUnit): 

556 vars.add(m) 

557 

558 # we only want radicals, so exclude Number handling; in this case 

559 # d will be evaluated 

560 d = collect_const(expr, *vars, Numbers=False) 

561 hit = expr != d 

562 

563 if not evaluate: 

564 nrad = 0 

565 # make the evaluated args canonical 

566 args = list(ordered(Add.make_args(d))) 

567 for i, m in enumerate(args): 

568 c, nc = m.args_cnc() 

569 for ci in c: 

570 # XXX should this be restricted to ci.is_number as above? 

571 if ci.is_Pow and ci.exp.is_Rational and ci.exp.q == 2 or \ 

572 ci is S.ImaginaryUnit: 

573 nrad += 1 

574 break 

575 args[i] *= coeff 

576 if not (hit or nrad): 

577 args = [Add(*args)] 

578 return tuple(args), nrad 

579 

580 return coeff*d 

581 

582 

583def collect_abs(expr): 

584 """Return ``expr`` with arguments of multiple Abs in a term collected 

585 under a single instance. 

586 

587 Examples 

588 ======== 

589 

590 >>> from sympy.simplify.radsimp import collect_abs 

591 >>> from sympy.abc import x 

592 >>> collect_abs(abs(x + 1)/abs(x**2 - 1)) 

593 Abs((x + 1)/(x**2 - 1)) 

594 >>> collect_abs(abs(1/x)) 

595 Abs(1/x) 

596 """ 

597 def _abs(mul): 

598 c, nc = mul.args_cnc() 

599 a = [] 

600 o = [] 

601 for i in c: 

602 if isinstance(i, Abs): 

603 a.append(i.args[0]) 

604 elif isinstance(i, Pow) and isinstance(i.base, Abs) and i.exp.is_real: 

605 a.append(i.base.args[0]**i.exp) 

606 else: 

607 o.append(i) 

608 if len(a) < 2 and not any(i.exp.is_negative for i in a if isinstance(i, Pow)): 

609 return mul 

610 absarg = Mul(*a) 

611 A = Abs(absarg) 

612 args = [A] 

613 args.extend(o) 

614 if not A.has(Abs): 

615 args.extend(nc) 

616 return Mul(*args) 

617 if not isinstance(A, Abs): 

618 # reevaluate and make it unevaluated 

619 A = Abs(absarg, evaluate=False) 

620 args[0] = A 

621 _mulsort(args) 

622 args.extend(nc) # nc always go last 

623 return Mul._from_args(args, is_commutative=not nc) 

624 

625 return expr.replace( 

626 lambda x: isinstance(x, Mul), 

627 lambda x: _abs(x)).replace( 

628 lambda x: isinstance(x, Pow), 

629 lambda x: _abs(x)) 

630 

631 

632def collect_const(expr, *vars, Numbers=True): 

633 """A non-greedy collection of terms with similar number coefficients in 

634 an Add expr. If ``vars`` is given then only those constants will be 

635 targeted. Although any Number can also be targeted, if this is not 

636 desired set ``Numbers=False`` and no Float or Rational will be collected. 

637 

638 Parameters 

639 ========== 

640 

641 expr : SymPy expression 

642 This parameter defines the expression the expression from which 

643 terms with similar coefficients are to be collected. A non-Add 

644 expression is returned as it is. 

645 

646 vars : variable length collection of Numbers, optional 

647 Specifies the constants to target for collection. Can be multiple in 

648 number. 

649 

650 Numbers : bool 

651 Specifies to target all instance of 

652 :class:`sympy.core.numbers.Number` class. If ``Numbers=False``, then 

653 no Float or Rational will be collected. 

654 

655 Returns 

656 ======= 

657 

658 expr : Expr 

659 Returns an expression with similar coefficient terms collected. 

660 

661 Examples 

662 ======== 

663 

664 >>> from sympy import sqrt 

665 >>> from sympy.abc import s, x, y, z 

666 >>> from sympy.simplify.radsimp import collect_const 

667 >>> collect_const(sqrt(3) + sqrt(3)*(1 + sqrt(2))) 

668 sqrt(3)*(sqrt(2) + 2) 

669 >>> collect_const(sqrt(3)*s + sqrt(7)*s + sqrt(3) + sqrt(7)) 

670 (sqrt(3) + sqrt(7))*(s + 1) 

671 >>> s = sqrt(2) + 2 

672 >>> collect_const(sqrt(3)*s + sqrt(3) + sqrt(7)*s + sqrt(7)) 

673 (sqrt(2) + 3)*(sqrt(3) + sqrt(7)) 

674 >>> collect_const(sqrt(3)*s + sqrt(3) + sqrt(7)*s + sqrt(7), sqrt(3)) 

675 sqrt(7) + sqrt(3)*(sqrt(2) + 3) + sqrt(7)*(sqrt(2) + 2) 

676 

677 The collection is sign-sensitive, giving higher precedence to the 

678 unsigned values: 

679 

680 >>> collect_const(x - y - z) 

681 x - (y + z) 

682 >>> collect_const(-y - z) 

683 -(y + z) 

684 >>> collect_const(2*x - 2*y - 2*z, 2) 

685 2*(x - y - z) 

686 >>> collect_const(2*x - 2*y - 2*z, -2) 

687 2*x - 2*(y + z) 

688 

689 See Also 

690 ======== 

691 

692 collect, collect_sqrt, rcollect 

693 """ 

694 if not expr.is_Add: 

695 return expr 

696 

697 recurse = False 

698 

699 if not vars: 

700 recurse = True 

701 vars = set() 

702 for a in expr.args: 

703 for m in Mul.make_args(a): 

704 if m.is_number: 

705 vars.add(m) 

706 else: 

707 vars = sympify(vars) 

708 if not Numbers: 

709 vars = [v for v in vars if not v.is_Number] 

710 

711 vars = list(ordered(vars)) 

712 for v in vars: 

713 terms = defaultdict(list) 

714 Fv = Factors(v) 

715 for m in Add.make_args(expr): 

716 f = Factors(m) 

717 q, r = f.div(Fv) 

718 if r.is_one: 

719 # only accept this as a true factor if 

720 # it didn't change an exponent from an Integer 

721 # to a non-Integer, e.g. 2/sqrt(2) -> sqrt(2) 

722 # -- we aren't looking for this sort of change 

723 fwas = f.factors.copy() 

724 fnow = q.factors 

725 if not any(k in fwas and fwas[k].is_Integer and not 

726 fnow[k].is_Integer for k in fnow): 

727 terms[v].append(q.as_expr()) 

728 continue 

729 terms[S.One].append(m) 

730 

731 args = [] 

732 hit = False 

733 uneval = False 

734 for k in ordered(terms): 

735 v = terms[k] 

736 if k is S.One: 

737 args.extend(v) 

738 continue 

739 

740 if len(v) > 1: 

741 v = Add(*v) 

742 hit = True 

743 if recurse and v != expr: 

744 vars.append(v) 

745 else: 

746 v = v[0] 

747 

748 # be careful not to let uneval become True unless 

749 # it must be because it's going to be more expensive 

750 # to rebuild the expression as an unevaluated one 

751 if Numbers and k.is_Number and v.is_Add: 

752 args.append(_keep_coeff(k, v, sign=True)) 

753 uneval = True 

754 else: 

755 args.append(k*v) 

756 

757 if hit: 

758 if uneval: 

759 expr = _unevaluated_Add(*args) 

760 else: 

761 expr = Add(*args) 

762 if not expr.is_Add: 

763 break 

764 

765 return expr 

766 

767 

768def radsimp(expr, symbolic=True, max_terms=4): 

769 r""" 

770 Rationalize the denominator by removing square roots. 

771 

772 Explanation 

773 =========== 

774 

775 The expression returned from radsimp must be used with caution 

776 since if the denominator contains symbols, it will be possible to make 

777 substitutions that violate the assumptions of the simplification process: 

778 that for a denominator matching a + b*sqrt(c), a != +/-b*sqrt(c). (If 

779 there are no symbols, this assumptions is made valid by collecting terms 

780 of sqrt(c) so the match variable ``a`` does not contain ``sqrt(c)``.) If 

781 you do not want the simplification to occur for symbolic denominators, set 

782 ``symbolic`` to False. 

783 

784 If there are more than ``max_terms`` radical terms then the expression is 

785 returned unchanged. 

786 

787 Examples 

788 ======== 

789 

790 >>> from sympy import radsimp, sqrt, Symbol, pprint 

791 >>> from sympy import factor_terms, fraction, signsimp 

792 >>> from sympy.simplify.radsimp import collect_sqrt 

793 >>> from sympy.abc import a, b, c 

794 

795 >>> radsimp(1/(2 + sqrt(2))) 

796 (2 - sqrt(2))/2 

797 >>> x,y = map(Symbol, 'xy') 

798 >>> e = ((2 + 2*sqrt(2))*x + (2 + sqrt(8))*y)/(2 + sqrt(2)) 

799 >>> radsimp(e) 

800 sqrt(2)*(x + y) 

801 

802 No simplification beyond removal of the gcd is done. One might 

803 want to polish the result a little, however, by collecting 

804 square root terms: 

805 

806 >>> r2 = sqrt(2) 

807 >>> r5 = sqrt(5) 

808 >>> ans = radsimp(1/(y*r2 + x*r2 + a*r5 + b*r5)); pprint(ans) 

809 ___ ___ ___ ___ 

810 \/ 5 *a + \/ 5 *b - \/ 2 *x - \/ 2 *y 

811 ------------------------------------------ 

812 2 2 2 2 

813 5*a + 10*a*b + 5*b - 2*x - 4*x*y - 2*y 

814 

815 >>> n, d = fraction(ans) 

816 >>> pprint(factor_terms(signsimp(collect_sqrt(n))/d, radical=True)) 

817 ___ ___ 

818 \/ 5 *(a + b) - \/ 2 *(x + y) 

819 ------------------------------------------ 

820 2 2 2 2 

821 5*a + 10*a*b + 5*b - 2*x - 4*x*y - 2*y 

822 

823 If radicals in the denominator cannot be removed or there is no denominator, 

824 the original expression will be returned. 

825 

826 >>> radsimp(sqrt(2)*x + sqrt(2)) 

827 sqrt(2)*x + sqrt(2) 

828 

829 Results with symbols will not always be valid for all substitutions: 

830 

831 >>> eq = 1/(a + b*sqrt(c)) 

832 >>> eq.subs(a, b*sqrt(c)) 

833 1/(2*b*sqrt(c)) 

834 >>> radsimp(eq).subs(a, b*sqrt(c)) 

835 nan 

836 

837 If ``symbolic=False``, symbolic denominators will not be transformed (but 

838 numeric denominators will still be processed): 

839 

840 >>> radsimp(eq, symbolic=False) 

841 1/(a + b*sqrt(c)) 

842 

843 """ 

844 from sympy.simplify.simplify import signsimp 

845 

846 syms = symbols("a:d A:D") 

847 def _num(rterms): 

848 # return the multiplier that will simplify the expression described 

849 # by rterms [(sqrt arg, coeff), ... ] 

850 a, b, c, d, A, B, C, D = syms 

851 if len(rterms) == 2: 

852 reps = dict(list(zip([A, a, B, b], [j for i in rterms for j in i]))) 

853 return ( 

854 sqrt(A)*a - sqrt(B)*b).xreplace(reps) 

855 if len(rterms) == 3: 

856 reps = dict(list(zip([A, a, B, b, C, c], [j for i in rterms for j in i]))) 

857 return ( 

858 (sqrt(A)*a + sqrt(B)*b - sqrt(C)*c)*(2*sqrt(A)*sqrt(B)*a*b - A*a**2 - 

859 B*b**2 + C*c**2)).xreplace(reps) 

860 elif len(rterms) == 4: 

861 reps = dict(list(zip([A, a, B, b, C, c, D, d], [j for i in rterms for j in i]))) 

862 return ((sqrt(A)*a + sqrt(B)*b - sqrt(C)*c - sqrt(D)*d)*(2*sqrt(A)*sqrt(B)*a*b 

863 - A*a**2 - B*b**2 - 2*sqrt(C)*sqrt(D)*c*d + C*c**2 + 

864 D*d**2)*(-8*sqrt(A)*sqrt(B)*sqrt(C)*sqrt(D)*a*b*c*d + A**2*a**4 - 

865 2*A*B*a**2*b**2 - 2*A*C*a**2*c**2 - 2*A*D*a**2*d**2 + B**2*b**4 - 

866 2*B*C*b**2*c**2 - 2*B*D*b**2*d**2 + C**2*c**4 - 2*C*D*c**2*d**2 + 

867 D**2*d**4)).xreplace(reps) 

868 elif len(rterms) == 1: 

869 return sqrt(rterms[0][0]) 

870 else: 

871 raise NotImplementedError 

872 

873 def ispow2(d, log2=False): 

874 if not d.is_Pow: 

875 return False 

876 e = d.exp 

877 if e.is_Rational and e.q == 2 or symbolic and denom(e) == 2: 

878 return True 

879 if log2: 

880 q = 1 

881 if e.is_Rational: 

882 q = e.q 

883 elif symbolic: 

884 d = denom(e) 

885 if d.is_Integer: 

886 q = d 

887 if q != 1 and log(q, 2).is_Integer: 

888 return True 

889 return False 

890 

891 def handle(expr): 

892 # Handle first reduces to the case 

893 # expr = 1/d, where d is an add, or d is base**p/2. 

894 # We do this by recursively calling handle on each piece. 

895 from sympy.simplify.simplify import nsimplify 

896 

897 n, d = fraction(expr) 

898 

899 if expr.is_Atom or (d.is_Atom and n.is_Atom): 

900 return expr 

901 elif not n.is_Atom: 

902 n = n.func(*[handle(a) for a in n.args]) 

903 return _unevaluated_Mul(n, handle(1/d)) 

904 elif n is not S.One: 

905 return _unevaluated_Mul(n, handle(1/d)) 

906 elif d.is_Mul: 

907 return _unevaluated_Mul(*[handle(1/d) for d in d.args]) 

908 

909 # By this step, expr is 1/d, and d is not a mul. 

910 if not symbolic and d.free_symbols: 

911 return expr 

912 

913 if ispow2(d): 

914 d2 = sqrtdenest(sqrt(d.base))**numer(d.exp) 

915 if d2 != d: 

916 return handle(1/d2) 

917 elif d.is_Pow and (d.exp.is_integer or d.base.is_positive): 

918 # (1/d**i) = (1/d)**i 

919 return handle(1/d.base)**d.exp 

920 

921 if not (d.is_Add or ispow2(d)): 

922 return 1/d.func(*[handle(a) for a in d.args]) 

923 

924 # handle 1/d treating d as an Add (though it may not be) 

925 

926 keep = True # keep changes that are made 

927 

928 # flatten it and collect radicals after checking for special 

929 # conditions 

930 d = _mexpand(d) 

931 

932 # did it change? 

933 if d.is_Atom: 

934 return 1/d 

935 

936 # is it a number that might be handled easily? 

937 if d.is_number: 

938 _d = nsimplify(d) 

939 if _d.is_Number and _d.equals(d): 

940 return 1/_d 

941 

942 while True: 

943 # collect similar terms 

944 collected = defaultdict(list) 

945 for m in Add.make_args(d): # d might have become non-Add 

946 p2 = [] 

947 other = [] 

948 for i in Mul.make_args(m): 

949 if ispow2(i, log2=True): 

950 p2.append(i.base if i.exp is S.Half else i.base**(2*i.exp)) 

951 elif i is S.ImaginaryUnit: 

952 p2.append(S.NegativeOne) 

953 else: 

954 other.append(i) 

955 collected[tuple(ordered(p2))].append(Mul(*other)) 

956 rterms = list(ordered(list(collected.items()))) 

957 rterms = [(Mul(*i), Add(*j)) for i, j in rterms] 

958 nrad = len(rterms) - (1 if rterms[0][0] is S.One else 0) 

959 if nrad < 1: 

960 break 

961 elif nrad > max_terms: 

962 # there may have been invalid operations leading to this point 

963 # so don't keep changes, e.g. this expression is troublesome 

964 # in collecting terms so as not to raise the issue of 2834: 

965 # r = sqrt(sqrt(5) + 5) 

966 # eq = 1/(sqrt(5)*r + 2*sqrt(5)*sqrt(-sqrt(5) + 5) + 5*r) 

967 keep = False 

968 break 

969 if len(rterms) > 4: 

970 # in general, only 4 terms can be removed with repeated squaring 

971 # but other considerations can guide selection of radical terms 

972 # so that radicals are removed 

973 if all(x.is_Integer and (y**2).is_Rational for x, y in rterms): 

974 nd, d = rad_rationalize(S.One, Add._from_args( 

975 [sqrt(x)*y for x, y in rterms])) 

976 n *= nd 

977 else: 

978 # is there anything else that might be attempted? 

979 keep = False 

980 break 

981 from sympy.simplify.powsimp import powsimp, powdenest 

982 

983 num = powsimp(_num(rterms)) 

984 n *= num 

985 d *= num 

986 d = powdenest(_mexpand(d), force=symbolic) 

987 if d.has(S.Zero, nan, zoo): 

988 return expr 

989 if d.is_Atom: 

990 break 

991 

992 if not keep: 

993 return expr 

994 return _unevaluated_Mul(n, 1/d) 

995 

996 coeff, expr = expr.as_coeff_Add() 

997 expr = expr.normal() 

998 old = fraction(expr) 

999 n, d = fraction(handle(expr)) 

1000 if old != (n, d): 

1001 if not d.is_Atom: 

1002 was = (n, d) 

1003 n = signsimp(n, evaluate=False) 

1004 d = signsimp(d, evaluate=False) 

1005 u = Factors(_unevaluated_Mul(n, 1/d)) 

1006 u = _unevaluated_Mul(*[k**v for k, v in u.factors.items()]) 

1007 n, d = fraction(u) 

1008 if old == (n, d): 

1009 n, d = was 

1010 n = expand_mul(n) 

1011 if d.is_Number or d.is_Add: 

1012 n2, d2 = fraction(gcd_terms(_unevaluated_Mul(n, 1/d))) 

1013 if d2.is_Number or (d2.count_ops() <= d.count_ops()): 

1014 n, d = [signsimp(i) for i in (n2, d2)] 

1015 if n.is_Mul and n.args[0].is_Number: 

1016 n = n.func(*n.args) 

1017 

1018 return coeff + _unevaluated_Mul(n, 1/d) 

1019 

1020 

1021def rad_rationalize(num, den): 

1022 """ 

1023 Rationalize ``num/den`` by removing square roots in the denominator; 

1024 num and den are sum of terms whose squares are positive rationals. 

1025 

1026 Examples 

1027 ======== 

1028 

1029 >>> from sympy import sqrt 

1030 >>> from sympy.simplify.radsimp import rad_rationalize 

1031 >>> rad_rationalize(sqrt(3), 1 + sqrt(2)/3) 

1032 (-sqrt(3) + sqrt(6)/3, -7/9) 

1033 """ 

1034 if not den.is_Add: 

1035 return num, den 

1036 g, a, b = split_surds(den) 

1037 a = a*sqrt(g) 

1038 num = _mexpand((a - b)*num) 

1039 den = _mexpand(a**2 - b**2) 

1040 return rad_rationalize(num, den) 

1041 

1042 

1043def fraction(expr, exact=False): 

1044 """Returns a pair with expression's numerator and denominator. 

1045 If the given expression is not a fraction then this function 

1046 will return the tuple (expr, 1). 

1047 

1048 This function will not make any attempt to simplify nested 

1049 fractions or to do any term rewriting at all. 

1050 

1051 If only one of the numerator/denominator pair is needed then 

1052 use numer(expr) or denom(expr) functions respectively. 

1053 

1054 >>> from sympy import fraction, Rational, Symbol 

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

1056 

1057 >>> fraction(x/y) 

1058 (x, y) 

1059 >>> fraction(x) 

1060 (x, 1) 

1061 

1062 >>> fraction(1/y**2) 

1063 (1, y**2) 

1064 

1065 >>> fraction(x*y/2) 

1066 (x*y, 2) 

1067 >>> fraction(Rational(1, 2)) 

1068 (1, 2) 

1069 

1070 This function will also work fine with assumptions: 

1071 

1072 >>> k = Symbol('k', negative=True) 

1073 >>> fraction(x * y**k) 

1074 (x, y**(-k)) 

1075 

1076 If we know nothing about sign of some exponent and ``exact`` 

1077 flag is unset, then structure this exponent's structure will 

1078 be analyzed and pretty fraction will be returned: 

1079 

1080 >>> from sympy import exp, Mul 

1081 >>> fraction(2*x**(-y)) 

1082 (2, x**y) 

1083 

1084 >>> fraction(exp(-x)) 

1085 (1, exp(x)) 

1086 

1087 >>> fraction(exp(-x), exact=True) 

1088 (exp(-x), 1) 

1089 

1090 The ``exact`` flag will also keep any unevaluated Muls from 

1091 being evaluated: 

1092 

1093 >>> u = Mul(2, x + 1, evaluate=False) 

1094 >>> fraction(u) 

1095 (2*x + 2, 1) 

1096 >>> fraction(u, exact=True) 

1097 (2*(x + 1), 1) 

1098 """ 

1099 expr = sympify(expr) 

1100 

1101 numer, denom = [], [] 

1102 

1103 for term in Mul.make_args(expr): 

1104 if term.is_commutative and (term.is_Pow or isinstance(term, exp)): 

1105 b, ex = term.as_base_exp() 

1106 if ex.is_negative: 

1107 if ex is S.NegativeOne: 

1108 denom.append(b) 

1109 elif exact: 

1110 if ex.is_constant(): 

1111 denom.append(Pow(b, -ex)) 

1112 else: 

1113 numer.append(term) 

1114 else: 

1115 denom.append(Pow(b, -ex)) 

1116 elif ex.is_positive: 

1117 numer.append(term) 

1118 elif not exact and ex.is_Mul: 

1119 n, d = term.as_numer_denom() 

1120 if n != 1: 

1121 numer.append(n) 

1122 denom.append(d) 

1123 else: 

1124 numer.append(term) 

1125 elif term.is_Rational and not term.is_Integer: 

1126 if term.p != 1: 

1127 numer.append(term.p) 

1128 denom.append(term.q) 

1129 else: 

1130 numer.append(term) 

1131 return Mul(*numer, evaluate=not exact), Mul(*denom, evaluate=not exact) 

1132 

1133 

1134def numer(expr): 

1135 return fraction(expr)[0] 

1136 

1137 

1138def denom(expr): 

1139 return fraction(expr)[1] 

1140 

1141 

1142def fraction_expand(expr, **hints): 

1143 return expr.expand(frac=True, **hints) 

1144 

1145 

1146def numer_expand(expr, **hints): 

1147 a, b = fraction(expr) 

1148 return a.expand(numer=True, **hints) / b 

1149 

1150 

1151def denom_expand(expr, **hints): 

1152 a, b = fraction(expr) 

1153 return a / b.expand(denom=True, **hints) 

1154 

1155 

1156expand_numer = numer_expand 

1157expand_denom = denom_expand 

1158expand_fraction = fraction_expand 

1159 

1160 

1161def split_surds(expr): 

1162 """ 

1163 Split an expression with terms whose squares are positive rationals 

1164 into a sum of terms whose surds squared have gcd equal to g 

1165 and a sum of terms with surds squared prime with g. 

1166 

1167 Examples 

1168 ======== 

1169 

1170 >>> from sympy import sqrt 

1171 >>> from sympy.simplify.radsimp import split_surds 

1172 >>> split_surds(3*sqrt(3) + sqrt(5)/7 + sqrt(6) + sqrt(10) + sqrt(15)) 

1173 (3, sqrt(2) + sqrt(5) + 3, sqrt(5)/7 + sqrt(10)) 

1174 """ 

1175 args = sorted(expr.args, key=default_sort_key) 

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

1177 surds = [x[1]**2 for x in coeff_muls if x[1].is_Pow] 

1178 surds.sort(key=default_sort_key) 

1179 g, b1, b2 = _split_gcd(*surds) 

1180 g2 = g 

1181 if not b2 and len(b1) >= 2: 

1182 b1n = [x/g for x in b1] 

1183 b1n = [x for x in b1n if x != 1] 

1184 # only a common factor has been factored; split again 

1185 g1, b1n, b2 = _split_gcd(*b1n) 

1186 g2 = g*g1 

1187 a1v, a2v = [], [] 

1188 for c, s in coeff_muls: 

1189 if s.is_Pow and s.exp == S.Half: 

1190 s1 = s.base 

1191 if s1 in b1: 

1192 a1v.append(c*sqrt(s1/g2)) 

1193 else: 

1194 a2v.append(c*s) 

1195 else: 

1196 a2v.append(c*s) 

1197 a = Add(*a1v) 

1198 b = Add(*a2v) 

1199 return g2, a, b 

1200 

1201 

1202def _split_gcd(*a): 

1203 """ 

1204 Split the list of integers ``a`` into a list of integers, ``a1`` having 

1205 ``g = gcd(a1)``, and a list ``a2`` whose elements are not divisible by 

1206 ``g``. Returns ``g, a1, a2``. 

1207 

1208 Examples 

1209 ======== 

1210 

1211 >>> from sympy.simplify.radsimp import _split_gcd 

1212 >>> _split_gcd(55, 35, 22, 14, 77, 10) 

1213 (5, [55, 35, 10], [22, 14, 77]) 

1214 """ 

1215 g = a[0] 

1216 b1 = [g] 

1217 b2 = [] 

1218 for x in a[1:]: 

1219 g1 = gcd(g, x) 

1220 if g1 == 1: 

1221 b2.append(x) 

1222 else: 

1223 g = g1 

1224 b1.append(x) 

1225 return g, b1, b2