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

348 statements  

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

1from collections import defaultdict 

2from functools import reduce 

3from math import prod 

4 

5from sympy.core.function import expand_log, count_ops, _coeff_isneg 

6from sympy.core import sympify, Basic, Dummy, S, Add, Mul, Pow, expand_mul, factor_terms 

7from sympy.core.sorting import ordered, default_sort_key 

8from sympy.core.numbers import Integer, Rational 

9from sympy.core.mul import _keep_coeff 

10from sympy.core.rules import Transform 

11from sympy.functions import exp_polar, exp, log, root, polarify, unpolarify 

12from sympy.matrices.expressions.matexpr import MatrixSymbol 

13from sympy.polys import lcm, gcd 

14from sympy.ntheory.factor_ import multiplicity 

15 

16 

17 

18def powsimp(expr, deep=False, combine='all', force=False, measure=count_ops): 

19 """ 

20 Reduce expression by combining powers with similar bases and exponents. 

21 

22 Explanation 

23 =========== 

24 

25 If ``deep`` is ``True`` then powsimp() will also simplify arguments of 

26 functions. By default ``deep`` is set to ``False``. 

27 

28 If ``force`` is ``True`` then bases will be combined without checking for 

29 assumptions, e.g. sqrt(x)*sqrt(y) -> sqrt(x*y) which is not true 

30 if x and y are both negative. 

31 

32 You can make powsimp() only combine bases or only combine exponents by 

33 changing combine='base' or combine='exp'. By default, combine='all', 

34 which does both. combine='base' will only combine:: 

35 

36 a a a 2x x 

37 x * y => (x*y) as well as things like 2 => 4 

38 

39 and combine='exp' will only combine 

40 :: 

41 

42 a b (a + b) 

43 x * x => x 

44 

45 combine='exp' will strictly only combine exponents in the way that used 

46 to be automatic. Also use deep=True if you need the old behavior. 

47 

48 When combine='all', 'exp' is evaluated first. Consider the first 

49 example below for when there could be an ambiguity relating to this. 

50 This is done so things like the second example can be completely 

51 combined. If you want 'base' combined first, do something like 

52 powsimp(powsimp(expr, combine='base'), combine='exp'). 

53 

54 Examples 

55 ======== 

56 

57 >>> from sympy import powsimp, exp, log, symbols 

58 >>> from sympy.abc import x, y, z, n 

59 >>> powsimp(x**y*x**z*y**z, combine='all') 

60 x**(y + z)*y**z 

61 >>> powsimp(x**y*x**z*y**z, combine='exp') 

62 x**(y + z)*y**z 

63 >>> powsimp(x**y*x**z*y**z, combine='base', force=True) 

64 x**y*(x*y)**z 

65 

66 >>> powsimp(x**z*x**y*n**z*n**y, combine='all', force=True) 

67 (n*x)**(y + z) 

68 >>> powsimp(x**z*x**y*n**z*n**y, combine='exp') 

69 n**(y + z)*x**(y + z) 

70 >>> powsimp(x**z*x**y*n**z*n**y, combine='base', force=True) 

71 (n*x)**y*(n*x)**z 

72 

73 >>> x, y = symbols('x y', positive=True) 

74 >>> powsimp(log(exp(x)*exp(y))) 

75 log(exp(x)*exp(y)) 

76 >>> powsimp(log(exp(x)*exp(y)), deep=True) 

77 x + y 

78 

79 Radicals with Mul bases will be combined if combine='exp' 

80 

81 >>> from sympy import sqrt 

82 >>> x, y = symbols('x y') 

83 

84 Two radicals are automatically joined through Mul: 

85 

86 >>> a=sqrt(x*sqrt(y)) 

87 >>> a*a**3 == a**4 

88 True 

89 

90 But if an integer power of that radical has been 

91 autoexpanded then Mul does not join the resulting factors: 

92 

93 >>> a**4 # auto expands to a Mul, no longer a Pow 

94 x**2*y 

95 >>> _*a # so Mul doesn't combine them 

96 x**2*y*sqrt(x*sqrt(y)) 

97 >>> powsimp(_) # but powsimp will 

98 (x*sqrt(y))**(5/2) 

99 >>> powsimp(x*y*a) # but won't when doing so would violate assumptions 

100 x*y*sqrt(x*sqrt(y)) 

101 

102 """ 

103 def recurse(arg, **kwargs): 

104 _deep = kwargs.get('deep', deep) 

105 _combine = kwargs.get('combine', combine) 

106 _force = kwargs.get('force', force) 

107 _measure = kwargs.get('measure', measure) 

108 return powsimp(arg, _deep, _combine, _force, _measure) 

109 

110 expr = sympify(expr) 

111 

112 if (not isinstance(expr, Basic) or isinstance(expr, MatrixSymbol) or ( 

113 expr.is_Atom or expr in (exp_polar(0), exp_polar(1)))): 

114 return expr 

115 

116 if deep or expr.is_Add or expr.is_Mul and _y not in expr.args: 

117 expr = expr.func(*[recurse(w) for w in expr.args]) 

118 

119 if expr.is_Pow: 

120 return recurse(expr*_y, deep=False)/_y 

121 

122 if not expr.is_Mul: 

123 return expr 

124 

125 # handle the Mul 

126 if combine in ('exp', 'all'): 

127 # Collect base/exp data, while maintaining order in the 

128 # non-commutative parts of the product 

129 c_powers = defaultdict(list) 

130 nc_part = [] 

131 newexpr = [] 

132 coeff = S.One 

133 for term in expr.args: 

134 if term.is_Rational: 

135 coeff *= term 

136 continue 

137 if term.is_Pow: 

138 term = _denest_pow(term) 

139 if term.is_commutative: 

140 b, e = term.as_base_exp() 

141 if deep: 

142 b, e = [recurse(i) for i in [b, e]] 

143 if b.is_Pow or isinstance(b, exp): 

144 # don't let smthg like sqrt(x**a) split into x**a, 1/2 

145 # or else it will be joined as x**(a/2) later 

146 b, e = b**e, S.One 

147 c_powers[b].append(e) 

148 else: 

149 # This is the logic that combines exponents for equal, 

150 # but non-commutative bases: A**x*A**y == A**(x+y). 

151 if nc_part: 

152 b1, e1 = nc_part[-1].as_base_exp() 

153 b2, e2 = term.as_base_exp() 

154 if (b1 == b2 and 

155 e1.is_commutative and e2.is_commutative): 

156 nc_part[-1] = Pow(b1, Add(e1, e2)) 

157 continue 

158 nc_part.append(term) 

159 

160 # add up exponents of common bases 

161 for b, e in ordered(iter(c_powers.items())): 

162 # allow 2**x/4 -> 2**(x - 2); don't do this when b and e are 

163 # Numbers since autoevaluation will undo it, e.g. 

164 # 2**(1/3)/4 -> 2**(1/3 - 2) -> 2**(1/3)/4 

165 if (b and b.is_Rational and not all(ei.is_Number for ei in e) and \ 

166 coeff is not S.One and 

167 b not in (S.One, S.NegativeOne)): 

168 m = multiplicity(abs(b), abs(coeff)) 

169 if m: 

170 e.append(m) 

171 coeff /= b**m 

172 c_powers[b] = Add(*e) 

173 if coeff is not S.One: 

174 if coeff in c_powers: 

175 c_powers[coeff] += S.One 

176 else: 

177 c_powers[coeff] = S.One 

178 

179 # convert to plain dictionary 

180 c_powers = dict(c_powers) 

181 

182 # check for base and inverted base pairs 

183 be = list(c_powers.items()) 

184 skip = set() # skip if we already saw them 

185 for b, e in be: 

186 if b in skip: 

187 continue 

188 bpos = b.is_positive or b.is_polar 

189 if bpos: 

190 binv = 1/b 

191 if b != binv and binv in c_powers: 

192 if b.as_numer_denom()[0] is S.One: 

193 c_powers.pop(b) 

194 c_powers[binv] -= e 

195 else: 

196 skip.add(binv) 

197 e = c_powers.pop(binv) 

198 c_powers[b] -= e 

199 

200 # check for base and negated base pairs 

201 be = list(c_powers.items()) 

202 _n = S.NegativeOne 

203 for b, e in be: 

204 if (b.is_Symbol or b.is_Add) and -b in c_powers and b in c_powers: 

205 if (b.is_positive is not None or e.is_integer): 

206 if e.is_integer or b.is_negative: 

207 c_powers[-b] += c_powers.pop(b) 

208 else: # (-b).is_positive so use its e 

209 e = c_powers.pop(-b) 

210 c_powers[b] += e 

211 if _n in c_powers: 

212 c_powers[_n] += e 

213 else: 

214 c_powers[_n] = e 

215 

216 # filter c_powers and convert to a list 

217 c_powers = [(b, e) for b, e in c_powers.items() if e] 

218 

219 # ============================================================== 

220 # check for Mul bases of Rational powers that can be combined with 

221 # separated bases, e.g. x*sqrt(x*y)*sqrt(x*sqrt(x*y)) -> 

222 # (x*sqrt(x*y))**(3/2) 

223 # ---------------- helper functions 

224 

225 def ratq(x): 

226 '''Return Rational part of x's exponent as it appears in the bkey. 

227 ''' 

228 return bkey(x)[0][1] 

229 

230 def bkey(b, e=None): 

231 '''Return (b**s, c.q), c.p where e -> c*s. If e is not given then 

232 it will be taken by using as_base_exp() on the input b. 

233 e.g. 

234 x**3/2 -> (x, 2), 3 

235 x**y -> (x**y, 1), 1 

236 x**(2*y/3) -> (x**y, 3), 2 

237 exp(x/2) -> (exp(a), 2), 1 

238 

239 ''' 

240 if e is not None: # coming from c_powers or from below 

241 if e.is_Integer: 

242 return (b, S.One), e 

243 elif e.is_Rational: 

244 return (b, Integer(e.q)), Integer(e.p) 

245 else: 

246 c, m = e.as_coeff_Mul(rational=True) 

247 if c is not S.One: 

248 if m.is_integer: 

249 return (b, Integer(c.q)), m*Integer(c.p) 

250 return (b**m, Integer(c.q)), Integer(c.p) 

251 else: 

252 return (b**e, S.One), S.One 

253 else: 

254 return bkey(*b.as_base_exp()) 

255 

256 def update(b): 

257 '''Decide what to do with base, b. If its exponent is now an 

258 integer multiple of the Rational denominator, then remove it 

259 and put the factors of its base in the common_b dictionary or 

260 update the existing bases if necessary. If it has been zeroed 

261 out, simply remove the base. 

262 ''' 

263 newe, r = divmod(common_b[b], b[1]) 

264 if not r: 

265 common_b.pop(b) 

266 if newe: 

267 for m in Mul.make_args(b[0]**newe): 

268 b, e = bkey(m) 

269 if b not in common_b: 

270 common_b[b] = 0 

271 common_b[b] += e 

272 if b[1] != 1: 

273 bases.append(b) 

274 # ---------------- end of helper functions 

275 

276 # assemble a dictionary of the factors having a Rational power 

277 common_b = {} 

278 done = [] 

279 bases = [] 

280 for b, e in c_powers: 

281 b, e = bkey(b, e) 

282 if b in common_b: 

283 common_b[b] = common_b[b] + e 

284 else: 

285 common_b[b] = e 

286 if b[1] != 1 and b[0].is_Mul: 

287 bases.append(b) 

288 bases.sort(key=default_sort_key) # this makes tie-breaking canonical 

289 bases.sort(key=measure, reverse=True) # handle longest first 

290 for base in bases: 

291 if base not in common_b: # it may have been removed already 

292 continue 

293 b, exponent = base 

294 last = False # True when no factor of base is a radical 

295 qlcm = 1 # the lcm of the radical denominators 

296 while True: 

297 bstart = b 

298 qstart = qlcm 

299 

300 bb = [] # list of factors 

301 ee = [] # (factor's expo. and it's current value in common_b) 

302 for bi in Mul.make_args(b): 

303 bib, bie = bkey(bi) 

304 if bib not in common_b or common_b[bib] < bie: 

305 ee = bb = [] # failed 

306 break 

307 ee.append([bie, common_b[bib]]) 

308 bb.append(bib) 

309 if ee: 

310 # find the number of integral extractions possible 

311 # e.g. [(1, 2), (2, 2)] -> min(2/1, 2/2) -> 1 

312 min1 = ee[0][1]//ee[0][0] 

313 for i in range(1, len(ee)): 

314 rat = ee[i][1]//ee[i][0] 

315 if rat < 1: 

316 break 

317 min1 = min(min1, rat) 

318 else: 

319 # update base factor counts 

320 # e.g. if ee = [(2, 5), (3, 6)] then min1 = 2 

321 # and the new base counts will be 5-2*2 and 6-2*3 

322 for i in range(len(bb)): 

323 common_b[bb[i]] -= min1*ee[i][0] 

324 update(bb[i]) 

325 # update the count of the base 

326 # e.g. x**2*y*sqrt(x*sqrt(y)) the count of x*sqrt(y) 

327 # will increase by 4 to give bkey (x*sqrt(y), 2, 5) 

328 common_b[base] += min1*qstart*exponent 

329 if (last # no more radicals in base 

330 or len(common_b) == 1 # nothing left to join with 

331 or all(k[1] == 1 for k in common_b) # no rad's in common_b 

332 ): 

333 break 

334 # see what we can exponentiate base by to remove any radicals 

335 # so we know what to search for 

336 # e.g. if base were x**(1/2)*y**(1/3) then we should 

337 # exponentiate by 6 and look for powers of x and y in the ratio 

338 # of 2 to 3 

339 qlcm = lcm([ratq(bi) for bi in Mul.make_args(bstart)]) 

340 if qlcm == 1: 

341 break # we are done 

342 b = bstart**qlcm 

343 qlcm *= qstart 

344 if all(ratq(bi) == 1 for bi in Mul.make_args(b)): 

345 last = True # we are going to be done after this next pass 

346 # this base no longer can find anything to join with and 

347 # since it was longer than any other we are done with it 

348 b, q = base 

349 done.append((b, common_b.pop(base)*Rational(1, q))) 

350 

351 # update c_powers and get ready to continue with powsimp 

352 c_powers = done 

353 # there may be terms still in common_b that were bases that were 

354 # identified as needing processing, so remove those, too 

355 for (b, q), e in common_b.items(): 

356 if (b.is_Pow or isinstance(b, exp)) and \ 

357 q is not S.One and not b.exp.is_Rational: 

358 b, be = b.as_base_exp() 

359 b = b**(be/q) 

360 else: 

361 b = root(b, q) 

362 c_powers.append((b, e)) 

363 check = len(c_powers) 

364 c_powers = dict(c_powers) 

365 assert len(c_powers) == check # there should have been no duplicates 

366 # ============================================================== 

367 

368 # rebuild the expression 

369 newexpr = expr.func(*(newexpr + [Pow(b, e) for b, e in c_powers.items()])) 

370 if combine == 'exp': 

371 return expr.func(newexpr, expr.func(*nc_part)) 

372 else: 

373 return recurse(expr.func(*nc_part), combine='base') * \ 

374 recurse(newexpr, combine='base') 

375 

376 elif combine == 'base': 

377 

378 # Build c_powers and nc_part. These must both be lists not 

379 # dicts because exp's are not combined. 

380 c_powers = [] 

381 nc_part = [] 

382 for term in expr.args: 

383 if term.is_commutative: 

384 c_powers.append(list(term.as_base_exp())) 

385 else: 

386 nc_part.append(term) 

387 

388 # Pull out numerical coefficients from exponent if assumptions allow 

389 # e.g., 2**(2*x) => 4**x 

390 for i in range(len(c_powers)): 

391 b, e = c_powers[i] 

392 if not (all(x.is_nonnegative for x in b.as_numer_denom()) or e.is_integer or force or b.is_polar): 

393 continue 

394 exp_c, exp_t = e.as_coeff_Mul(rational=True) 

395 if exp_c is not S.One and exp_t is not S.One: 

396 c_powers[i] = [Pow(b, exp_c), exp_t] 

397 

398 # Combine bases whenever they have the same exponent and 

399 # assumptions allow 

400 # first gather the potential bases under the common exponent 

401 c_exp = defaultdict(list) 

402 for b, e in c_powers: 

403 if deep: 

404 e = recurse(e) 

405 if e.is_Add and (b.is_positive or e.is_integer): 

406 e = factor_terms(e) 

407 if _coeff_isneg(e): 

408 e = -e 

409 b = 1/b 

410 c_exp[e].append(b) 

411 del c_powers 

412 

413 # Merge back in the results of the above to form a new product 

414 c_powers = defaultdict(list) 

415 for e in c_exp: 

416 bases = c_exp[e] 

417 

418 # calculate the new base for e 

419 

420 if len(bases) == 1: 

421 new_base = bases[0] 

422 elif e.is_integer or force: 

423 new_base = expr.func(*bases) 

424 else: 

425 # see which ones can be joined 

426 unk = [] 

427 nonneg = [] 

428 neg = [] 

429 for bi in bases: 

430 if bi.is_negative: 

431 neg.append(bi) 

432 elif bi.is_nonnegative: 

433 nonneg.append(bi) 

434 elif bi.is_polar: 

435 nonneg.append( 

436 bi) # polar can be treated like non-negative 

437 else: 

438 unk.append(bi) 

439 if len(unk) == 1 and not neg or len(neg) == 1 and not unk: 

440 # a single neg or a single unk can join the rest 

441 nonneg.extend(unk + neg) 

442 unk = neg = [] 

443 elif neg: 

444 # their negative signs cancel in groups of 2*q if we know 

445 # that e = p/q else we have to treat them as unknown 

446 israt = False 

447 if e.is_Rational: 

448 israt = True 

449 else: 

450 p, d = e.as_numer_denom() 

451 if p.is_integer and d.is_integer: 

452 israt = True 

453 if israt: 

454 neg = [-w for w in neg] 

455 unk.extend([S.NegativeOne]*len(neg)) 

456 else: 

457 unk.extend(neg) 

458 neg = [] 

459 del israt 

460 

461 # these shouldn't be joined 

462 for b in unk: 

463 c_powers[b].append(e) 

464 # here is a new joined base 

465 new_base = expr.func(*(nonneg + neg)) 

466 # if there are positive parts they will just get separated 

467 # again unless some change is made 

468 

469 def _terms(e): 

470 # return the number of terms of this expression 

471 # when multiplied out -- assuming no joining of terms 

472 if e.is_Add: 

473 return sum([_terms(ai) for ai in e.args]) 

474 if e.is_Mul: 

475 return prod([_terms(mi) for mi in e.args]) 

476 return 1 

477 xnew_base = expand_mul(new_base, deep=False) 

478 if len(Add.make_args(xnew_base)) < _terms(new_base): 

479 new_base = factor_terms(xnew_base) 

480 

481 c_powers[new_base].append(e) 

482 

483 # break out the powers from c_powers now 

484 c_part = [Pow(b, ei) for b, e in c_powers.items() for ei in e] 

485 

486 # we're done 

487 return expr.func(*(c_part + nc_part)) 

488 

489 else: 

490 raise ValueError("combine must be one of ('all', 'exp', 'base').") 

491 

492 

493def powdenest(eq, force=False, polar=False): 

494 r""" 

495 Collect exponents on powers as assumptions allow. 

496 

497 Explanation 

498 =========== 

499 

500 Given ``(bb**be)**e``, this can be simplified as follows: 

501 * if ``bb`` is positive, or 

502 * ``e`` is an integer, or 

503 * ``|be| < 1`` then this simplifies to ``bb**(be*e)`` 

504 

505 Given a product of powers raised to a power, ``(bb1**be1 * 

506 bb2**be2...)**e``, simplification can be done as follows: 

507 

508 - if e is positive, the gcd of all bei can be joined with e; 

509 - all non-negative bb can be separated from those that are negative 

510 and their gcd can be joined with e; autosimplification already 

511 handles this separation. 

512 - integer factors from powers that have integers in the denominator 

513 of the exponent can be removed from any term and the gcd of such 

514 integers can be joined with e 

515 

516 Setting ``force`` to ``True`` will make symbols that are not explicitly 

517 negative behave as though they are positive, resulting in more 

518 denesting. 

519 

520 Setting ``polar`` to ``True`` will do simplifications on the Riemann surface of 

521 the logarithm, also resulting in more denestings. 

522 

523 When there are sums of logs in exp() then a product of powers may be 

524 obtained e.g. ``exp(3*(log(a) + 2*log(b)))`` - > ``a**3*b**6``. 

525 

526 Examples 

527 ======== 

528 

529 >>> from sympy.abc import a, b, x, y, z 

530 >>> from sympy import Symbol, exp, log, sqrt, symbols, powdenest 

531 

532 >>> powdenest((x**(2*a/3))**(3*x)) 

533 (x**(2*a/3))**(3*x) 

534 >>> powdenest(exp(3*x*log(2))) 

535 2**(3*x) 

536 

537 Assumptions may prevent expansion: 

538 

539 >>> powdenest(sqrt(x**2)) 

540 sqrt(x**2) 

541 

542 >>> p = symbols('p', positive=True) 

543 >>> powdenest(sqrt(p**2)) 

544 p 

545 

546 No other expansion is done. 

547 

548 >>> i, j = symbols('i,j', integer=True) 

549 >>> powdenest((x**x)**(i + j)) # -X-> (x**x)**i*(x**x)**j 

550 x**(x*(i + j)) 

551 

552 But exp() will be denested by moving all non-log terms outside of 

553 the function; this may result in the collapsing of the exp to a power 

554 with a different base: 

555 

556 >>> powdenest(exp(3*y*log(x))) 

557 x**(3*y) 

558 >>> powdenest(exp(y*(log(a) + log(b)))) 

559 (a*b)**y 

560 >>> powdenest(exp(3*(log(a) + log(b)))) 

561 a**3*b**3 

562 

563 If assumptions allow, symbols can also be moved to the outermost exponent: 

564 

565 >>> i = Symbol('i', integer=True) 

566 >>> powdenest(((x**(2*i))**(3*y))**x) 

567 ((x**(2*i))**(3*y))**x 

568 >>> powdenest(((x**(2*i))**(3*y))**x, force=True) 

569 x**(6*i*x*y) 

570 

571 >>> powdenest(((x**(2*a/3))**(3*y/i))**x) 

572 ((x**(2*a/3))**(3*y/i))**x 

573 >>> powdenest((x**(2*i)*y**(4*i))**z, force=True) 

574 (x*y**2)**(2*i*z) 

575 

576 >>> n = Symbol('n', negative=True) 

577 

578 >>> powdenest((x**i)**y, force=True) 

579 x**(i*y) 

580 >>> powdenest((n**i)**x, force=True) 

581 (n**i)**x 

582 

583 """ 

584 from sympy.simplify.simplify import posify 

585 

586 if force: 

587 def _denest(b, e): 

588 if not isinstance(b, (Pow, exp)): 

589 return b.is_positive, Pow(b, e, evaluate=False) 

590 return _denest(b.base, b.exp*e) 

591 reps = [] 

592 for p in eq.atoms(Pow, exp): 

593 if isinstance(p.base, (Pow, exp)): 

594 ok, dp = _denest(*p.args) 

595 if ok is not False: 

596 reps.append((p, dp)) 

597 if reps: 

598 eq = eq.subs(reps) 

599 eq, reps = posify(eq) 

600 return powdenest(eq, force=False, polar=polar).xreplace(reps) 

601 

602 if polar: 

603 eq, rep = polarify(eq) 

604 return unpolarify(powdenest(unpolarify(eq, exponents_only=True)), rep) 

605 

606 new = powsimp(eq) 

607 return new.xreplace(Transform( 

608 _denest_pow, filter=lambda m: m.is_Pow or isinstance(m, exp))) 

609 

610_y = Dummy('y') 

611 

612 

613def _denest_pow(eq): 

614 """ 

615 Denest powers. 

616 

617 This is a helper function for powdenest that performs the actual 

618 transformation. 

619 """ 

620 from sympy.simplify.simplify import logcombine 

621 

622 b, e = eq.as_base_exp() 

623 if b.is_Pow or isinstance(b, exp) and e != 1: 

624 new = b._eval_power(e) 

625 if new is not None: 

626 eq = new 

627 b, e = new.as_base_exp() 

628 

629 # denest exp with log terms in exponent 

630 if b is S.Exp1 and e.is_Mul: 

631 logs = [] 

632 other = [] 

633 for ei in e.args: 

634 if any(isinstance(ai, log) for ai in Add.make_args(ei)): 

635 logs.append(ei) 

636 else: 

637 other.append(ei) 

638 logs = logcombine(Mul(*logs)) 

639 return Pow(exp(logs), Mul(*other)) 

640 

641 _, be = b.as_base_exp() 

642 if be is S.One and not (b.is_Mul or 

643 b.is_Rational and b.q != 1 or 

644 b.is_positive): 

645 return eq 

646 

647 # denest eq which is either pos**e or Pow**e or Mul**e or 

648 # Mul(b1**e1, b2**e2) 

649 

650 # handle polar numbers specially 

651 polars, nonpolars = [], [] 

652 for bb in Mul.make_args(b): 

653 if bb.is_polar: 

654 polars.append(bb.as_base_exp()) 

655 else: 

656 nonpolars.append(bb) 

657 if len(polars) == 1 and not polars[0][0].is_Mul: 

658 return Pow(polars[0][0], polars[0][1]*e)*powdenest(Mul(*nonpolars)**e) 

659 elif polars: 

660 return Mul(*[powdenest(bb**(ee*e)) for (bb, ee) in polars]) \ 

661 *powdenest(Mul(*nonpolars)**e) 

662 

663 if b.is_Integer: 

664 # use log to see if there is a power here 

665 logb = expand_log(log(b)) 

666 if logb.is_Mul: 

667 c, logb = logb.args 

668 e *= c 

669 base = logb.args[0] 

670 return Pow(base, e) 

671 

672 # if b is not a Mul or any factor is an atom then there is nothing to do 

673 if not b.is_Mul or any(s.is_Atom for s in Mul.make_args(b)): 

674 return eq 

675 

676 # let log handle the case of the base of the argument being a Mul, e.g. 

677 # sqrt(x**(2*i)*y**(6*i)) -> x**i*y**(3**i) if x and y are positive; we 

678 # will take the log, expand it, and then factor out the common powers that 

679 # now appear as coefficient. We do this manually since terms_gcd pulls out 

680 # fractions, terms_gcd(x+x*y/2) -> x*(y + 2)/2 and we don't want the 1/2; 

681 # gcd won't pull out numerators from a fraction: gcd(3*x, 9*x/2) -> x but 

682 # we want 3*x. Neither work with noncommutatives. 

683 

684 def nc_gcd(aa, bb): 

685 a, b = [i.as_coeff_Mul() for i in [aa, bb]] 

686 c = gcd(a[0], b[0]).as_numer_denom()[0] 

687 g = Mul(*(a[1].args_cnc(cset=True)[0] & b[1].args_cnc(cset=True)[0])) 

688 return _keep_coeff(c, g) 

689 

690 glogb = expand_log(log(b)) 

691 if glogb.is_Add: 

692 args = glogb.args 

693 g = reduce(nc_gcd, args) 

694 if g != 1: 

695 cg, rg = g.as_coeff_Mul() 

696 glogb = _keep_coeff(cg, rg*Add(*[a/g for a in args])) 

697 

698 # now put the log back together again 

699 if isinstance(glogb, log) or not glogb.is_Mul: 

700 if glogb.args[0].is_Pow or isinstance(glogb.args[0], exp): 

701 glogb = _denest_pow(glogb.args[0]) 

702 if (abs(glogb.exp) < 1) == True: 

703 return Pow(glogb.base, glogb.exp*e) 

704 return eq 

705 

706 # the log(b) was a Mul so join any adds with logcombine 

707 add = [] 

708 other = [] 

709 for a in glogb.args: 

710 if a.is_Add: 

711 add.append(a) 

712 else: 

713 other.append(a) 

714 return Pow(exp(logcombine(Mul(*add))), e*Mul(*other))