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
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from collections import defaultdict
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
22def collect(expr, syms, func=None, evaluate=None, exact=False, distribute_order_term=True):
23 """
24 Collect additive terms of an expression.
26 Explanation
27 ===========
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.
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.
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.
47 Examples
48 ========
50 >>> from sympy import S, collect, expand, factor, Wild
51 >>> from sympy.abc import a, b, c, x, y
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::
57 >>> collect(a*x**2 + b*x**2 + a*x - b*x + c, x)
58 c + x**2*(a + b) + x*(a - b)
60 The same result can be achieved in dictionary form::
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
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::
74 >>> collect(x**2 + y*x**2 + x*y + y + a*y, [x, y])
75 x**2*(y + 1) + x*y + y*(a + 1)
77 Also more complicated expressions can be used as patterns::
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)
83 >>> collect(a*x*log(x) + b*(x*log(x)), x*log(x))
84 x*(a + b)*log(x)
86 You can use wildcards in the pattern::
88 >>> w = Wild('w1')
89 >>> collect(a*x**y - b*x**y, w**y)
90 x**y*(a - b)
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::
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)
101 However if you incorporate rationals to the exponents, then you will get
102 well known behavior::
104 >>> collect(a*x**(2*c) + b*x**(2*c), x**c)
105 x**(2*c)*(a + b)
107 Note also that all previously stated facts about :func:`collect` function
108 apply to the exponential function, so you can get::
110 >>> from sympy import exp
111 >>> collect(a*exp(2*x) + b*exp(2*x), exp(x))
112 (a + b)*exp(2*x)
114 If you are interested only in collecting specific powers of some symbols
115 then set ``exact`` flag to True::
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)
122 If you want to collect on any object containing symbols, set
123 ``exact`` to None:
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)
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::
136 >>> from sympy import Derivative as D, collect, Function
137 >>> f = Function('f') (x)
139 >>> collect(a*D(f,x) + b*D(f,x), D(f,x))
140 (a + b)*Derivative(f(x), x)
142 >>> collect(a*D(D(f,x),x) + b*D(D(f,x),x), f)
143 (a + b)*Derivative(f(x), (x, 2))
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))
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)
151 Or you can even match both derivative order and exponent at the same time::
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
156 Finally, you can apply a function to each of the collected coefficients.
157 For example you can factorize symbolic coefficients of polynomial::
159 >>> f = expand((x + a + 1)**3)
161 >>> collect(f, x, factor)
162 x**3 + 3*x**2*(a + 1) + 3*x*(a + 1)**2 + (a + 1)**3
164 .. note:: Arguments are expected to be in expanded form, so you might have
165 to call :func:`~.expand` prior to calling this function.
167 See Also
168 ========
170 collect_const, collect_sqrt, rcollect
171 """
172 expr = sympify(expr)
173 syms = [sympify(i) for i in (syms if iterable(syms) else [syms])]
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()}
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)
212 if evaluate is None:
213 evaluate = global_parameters.evaluate
215 def make_expression(terms):
216 product = []
218 for term, rat, sym, deriv in terms:
219 if deriv is not None:
220 var, order = deriv
222 while order > 0:
223 term, order = Derivative(term, var), order - 1
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))
233 return Mul(*product)
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
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')
247 while isinstance(expr, Derivative):
248 s0 = expr.variables[0]
250 for s in expr.variables:
251 if s != s0:
252 raise NotImplementedError(
253 'Improve MV Derivative support in collect')
255 if s0 == sym:
256 expr, order = expr.expr, order + len(expr.variables)
257 else:
258 break
260 return expr, (sym, Rational(order))
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
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
277 if expr.is_Pow:
278 if isinstance(expr.base, Derivative):
279 sexpr, deriv = parse_derivative(expr.base)
280 else:
281 sexpr = expr.base
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
291 elif expr.exp.is_Number:
292 rat_expo = expr.exp
293 else:
294 coeff, tail = expr.exp.as_coeff_Mul()
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)
310 return sexpr, rat_expo, sym_expo, deriv
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)
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]
326 terms = terms[:] # need a copy
327 elems, common_expo, has_deriv = [], None, False
329 for elem, e_rat, e_sym, e_ord in pattern:
331 if elem.is_Number and e_rat == 1 and e_sym is None:
332 # a constant is a match for everything
333 continue
335 for j in range(len(terms)):
336 if terms[j] is None:
337 continue
339 term, t_rat, t_sym, t_ord = terms[j]
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
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
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
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
376 break
378 else:
379 # pattern element not found
380 return None
382 return [_f for _f in terms if _f], elems, common_expo, has_deriv
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)
399 syms = [expand_power_base(i, deep=False) for i in syms]
401 order_term = None
403 if distribute_order_term:
404 order_term = expr.getO()
406 if order_term is not None:
407 if order_term.has(*syms):
408 order_term = None
409 else:
410 expr = expr.removeO()
412 summa = [expand_power_base(i, deep=False) for i in Add.make_args(expr)]
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
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)
427 if result is not None:
428 if not symbol.is_commutative:
429 raise AttributeError("Can not collect noncommutative symbol")
431 terms, elems, common_expo, has_deriv = result
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()}
456 if disliked is not S.Zero:
457 collected[S.One] = disliked
459 if order_term is not None:
460 for key, val in collected.items():
461 collected[key] = val + order_term
463 if func is not None:
464 collected = {
465 key: func(val) for key, val in collected.items()}
467 if evaluate:
468 return Add(*[key*val for key, val in collected.items()])
469 else:
470 return collected
473def rcollect(expr, *vars):
474 """
475 Recursively collect sums in an expression.
477 Examples
478 ========
480 >>> from sympy.simplify import rcollect
481 >>> from sympy.abc import x, y
483 >>> expr = (x**2*y + x*y + x + y)/(x + y)
485 >>> rcollect(expr, y)
486 (x + y*(x**2 + x + 1))/(x + y)
488 See Also
489 ========
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])
498 if expr.is_Add:
499 return collect(expr, vars)
500 else:
501 return expr
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.
512 Note: since I = sqrt(-1), it is collected, too.
514 Examples
515 ========
517 >>> from sympy import sqrt
518 >>> from sympy.simplify.radsimp import collect_sqrt
519 >>> from sympy.abc import a, b
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)
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:
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)
540 See Also
541 ========
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)
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
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
580 return coeff*d
583def collect_abs(expr):
584 """Return ``expr`` with arguments of multiple Abs in a term collected
585 under a single instance.
587 Examples
588 ========
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)
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))
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.
638 Parameters
639 ==========
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.
646 vars : variable length collection of Numbers, optional
647 Specifies the constants to target for collection. Can be multiple in
648 number.
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.
655 Returns
656 =======
658 expr : Expr
659 Returns an expression with similar coefficient terms collected.
661 Examples
662 ========
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)
677 The collection is sign-sensitive, giving higher precedence to the
678 unsigned values:
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)
689 See Also
690 ========
692 collect, collect_sqrt, rcollect
693 """
694 if not expr.is_Add:
695 return expr
697 recurse = False
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]
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)
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
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]
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)
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
765 return expr
768def radsimp(expr, symbolic=True, max_terms=4):
769 r"""
770 Rationalize the denominator by removing square roots.
772 Explanation
773 ===========
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.
784 If there are more than ``max_terms`` radical terms then the expression is
785 returned unchanged.
787 Examples
788 ========
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
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)
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:
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
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
823 If radicals in the denominator cannot be removed or there is no denominator,
824 the original expression will be returned.
826 >>> radsimp(sqrt(2)*x + sqrt(2))
827 sqrt(2)*x + sqrt(2)
829 Results with symbols will not always be valid for all substitutions:
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
837 If ``symbolic=False``, symbolic denominators will not be transformed (but
838 numeric denominators will still be processed):
840 >>> radsimp(eq, symbolic=False)
841 1/(a + b*sqrt(c))
843 """
844 from sympy.simplify.simplify import signsimp
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
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
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
897 n, d = fraction(expr)
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])
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
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
921 if not (d.is_Add or ispow2(d)):
922 return 1/d.func(*[handle(a) for a in d.args])
924 # handle 1/d treating d as an Add (though it may not be)
926 keep = True # keep changes that are made
928 # flatten it and collect radicals after checking for special
929 # conditions
930 d = _mexpand(d)
932 # did it change?
933 if d.is_Atom:
934 return 1/d
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
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
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
992 if not keep:
993 return expr
994 return _unevaluated_Mul(n, 1/d)
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)
1018 return coeff + _unevaluated_Mul(n, 1/d)
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.
1026 Examples
1027 ========
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)
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).
1048 This function will not make any attempt to simplify nested
1049 fractions or to do any term rewriting at all.
1051 If only one of the numerator/denominator pair is needed then
1052 use numer(expr) or denom(expr) functions respectively.
1054 >>> from sympy import fraction, Rational, Symbol
1055 >>> from sympy.abc import x, y
1057 >>> fraction(x/y)
1058 (x, y)
1059 >>> fraction(x)
1060 (x, 1)
1062 >>> fraction(1/y**2)
1063 (1, y**2)
1065 >>> fraction(x*y/2)
1066 (x*y, 2)
1067 >>> fraction(Rational(1, 2))
1068 (1, 2)
1070 This function will also work fine with assumptions:
1072 >>> k = Symbol('k', negative=True)
1073 >>> fraction(x * y**k)
1074 (x, y**(-k))
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:
1080 >>> from sympy import exp, Mul
1081 >>> fraction(2*x**(-y))
1082 (2, x**y)
1084 >>> fraction(exp(-x))
1085 (1, exp(x))
1087 >>> fraction(exp(-x), exact=True)
1088 (exp(-x), 1)
1090 The ``exact`` flag will also keep any unevaluated Muls from
1091 being evaluated:
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)
1101 numer, denom = [], []
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)
1134def numer(expr):
1135 return fraction(expr)[0]
1138def denom(expr):
1139 return fraction(expr)[1]
1142def fraction_expand(expr, **hints):
1143 return expr.expand(frac=True, **hints)
1146def numer_expand(expr, **hints):
1147 a, b = fraction(expr)
1148 return a.expand(numer=True, **hints) / b
1151def denom_expand(expr, **hints):
1152 a, b = fraction(expr)
1153 return a / b.expand(denom=True, **hints)
1156expand_numer = numer_expand
1157expand_denom = denom_expand
1158expand_fraction = fraction_expand
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.
1167 Examples
1168 ========
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
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``.
1208 Examples
1209 ========
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