# Source code for sympy.core.mul

from basic import Basic, C
from singleton import S
from operations import AssocOp
from cache import cacheit
from logic import fuzzy_not
from compatibility import cmp_to_key

# internal marker to indicate:
#   "there are still non-commutative objects -- don't forget to process them"
class NC_Marker:
is_Order    = False
is_Mul      = False
is_Number   = False
is_Poly     = False

is_commutative = False

[docs]class Mul(AssocOp):

__slots__ = []

is_Mul = True

#identity = S.One
# cyclic import, so defined in numbers.py

@classmethod
def flatten(cls, seq):

# apply associativity, separate commutative part of seq
c_part = []         # out: commutative factors
nc_part = []        # out: non-commutative factors

nc_seq = []

coeff = S.One       # standalone term
# e.g. 3 * ...

c_powers = []       # (base,exp)      n
# e.g. (x,n) for x

num_exp = []        # (num-base, exp)           y
# e.g.  (3, y)  for  ... * 3  * ...

neg1e = 0           # exponent on -1 extracted from Number-based Pow

pnum_rat = {}       # (num-base, Rat-exp)          1/2
# e.g.  (3, 1/2)  for  ... * 3     * ...

order_symbols = None

# --- PART 1 ---
#
# "collect powers and coeff":
#
# o coeff
# o c_powers
# o num_exp
# o neg1e
# o pnum_rat
#
# NOTE: this is optimized for all-objects-are-commutative case

for o in seq:
# O(x)
if o.is_Order:
o, order_symbols = o.as_expr_variables(order_symbols)

# Mul([...])
if o.is_Mul:
if o.is_commutative:
seq.extend(o.args)    # XXX zerocopy?

else:
# NCMul can have commutative parts as well
for q in o.args:
if q.is_commutative:
seq.append(q)
else:
nc_seq.append(q)

# append non-commutative marker, so we don't forget to
# process scheduled non-commutative objects
seq.append(NC_Marker)

continue

# 3
elif o.is_Number:
if o is S.NaN or coeff is S.ComplexInfinity and o is S.Zero:
# we know for sure the result will be nan
return [S.NaN], [], None
elif coeff.is_Number: # it could be zoo
coeff *= o
if coeff is S.NaN:
# we know for sure the result will be nan
return [S.NaN], [], None
continue

elif o is S.ComplexInfinity:
if not coeff or coeff is S.ComplexInfinity:
# we know for sure the result will be nan
return [S.NaN], [], None
coeff = S.ComplexInfinity
continue

elif o.is_commutative:
#      e
# o = b
b, e = o.as_base_exp()

#  y
# 3
if o.is_Pow and b.is_Number:

# get all the factors with numeric base so they can be
# combined below, but don't combine negatives unless
# the exponent is an integer
if e.is_Rational:
if e.is_Integer:
coeff *= Pow(b, e) # it is an unevaluated power
continue
elif e.is_negative:    # also a sign of an unevaluated power
seq.append(Pow(b, e))
continue
elif b.is_negative:
neg1e += e
b = -b
if b is not S.One:
pnum_rat.setdefault(b, []).append(e)
continue
elif b.is_positive or e.is_integer:
num_exp.append((b, e))
continue
c_powers.append((b,e))

# NON-COMMUTATIVE
# TODO: Make non-commutative exponents not combine automatically
else:
if o is not NC_Marker:
nc_seq.append(o)

# process nc_seq (if any)
while nc_seq:
o = nc_seq.pop(0)
if not nc_part:
nc_part.append(o)
continue

#                             b    c       b+c
# try to combine last terms: a  * a   ->  a
o1 = nc_part.pop()
b1,e1 = o1.as_base_exp()
b2,e2 = o.as_base_exp()
new_exp = e1 + e2
# Only allow powers to combine if the new exponent is
# not an Add. This allow things like a**2*b**3 == a**5
# if a.is_commutative == False, but prohibits
# a**x*a**y and x**a*x**b from combining (x,y commute).
o12 = b1 ** new_exp

# now o12 could be a commutative object
if o12.is_commutative:
seq.append(o12)
continue
else:
nc_seq.insert(0, o12)

else:
nc_part.append(o1)
nc_part.append(o)

# We do want a combined exponent if it would not be an Add, such as
#  y    2y     3y
# x  * x   -> x
# We determine this if two exponents have the same term in as_coeff_mul
#
# Unfortunately, this isn't smart enough to consider combining into
#  z - y    y
# x      * x  will be left alone.  This is because checking every possible
# combination can slow things down.

# gather exponents of common bases...
# in c_powers
new_c_powers = []
common_b = {} # b:e
for b, e in c_powers:
co = e.as_coeff_mul()
common_b.setdefault(b, {}).setdefault(co[1], []).append(co[0])
for b, d in common_b.items():
for di, li in d.items():
for b, e in common_b.items():
for t, c in e.items():
new_c_powers.append((b,c*Mul(*t)))
c_powers = new_c_powers

# and in num_exp
new_num_exp = []
common_b = {} # b:e
for b, e in num_exp:
co = e.as_coeff_mul()
common_b.setdefault(b, {}).setdefault(co[1], []).append(co[0])
for b, d in common_b.items():
for di, li in d.items():
for b, e in common_b.items():
for t, c in e.items():
new_num_exp.append((b,c*Mul(*t)))
num_exp = new_num_exp

# --- PART 2 ---
#
# o process collected powers  (x**0 -> 1; x**1 -> x; otherwise Pow)
# o combine collected powers  (2**x * 3**x -> 6**x)
#   with numeric base

# ................................
# now we have:
# - coeff:
# - c_powers:    (b, e)
# - num_exp:     (2, e)
# - pnum_rat:    {(1/3, [1/3, 2/3, 1/4])}

#  0             1
# x  -> 1       x  -> x
for b, e in c_powers:
if e is S.One:
if b.is_Number:
coeff *= b
else:
c_part.append(b)
elif not e is S.Zero:
c_part.append(Pow(b, e))

#  x    x     x
# 2  * 3  -> 6
inv_exp_dict = {}   # exp:Mul(num-bases)     x    x
# e.g.  x:6  for  ... * 2  * 3  * ...
for b, e in num_exp:
inv_exp_dict.setdefault(e, []).append(b)
for e, b in inv_exp_dict.items():
inv_exp_dict[e] = Mul(*b)
c_part.extend([Pow(b, e) for e, b in inv_exp_dict.iteritems() if e])

# b, e -> e, b
# {(1/5, [1/3]), (1/2, [1/12, 1/4]} -> {(1/3, [1/5, 1/2])}
comb_e = {}
for b, e in pnum_rat.iteritems():
del pnum_rat
# process them, reducing exponents to values less than 1
# and updating coeff if necessary else adding them to
# num_rat for further processing
num_rat = []
for e, b in comb_e.iteritems():
b = Mul(*b)
if e.q == 1:
coeff *= Pow(b, e)
continue
if e.p > e.q:
e_i, ep = divmod(e.p, e.q)
coeff *= Pow(b, e_i)
e = Rational(ep, e.q)
num_rat.append((b, e))
del comb_e

# extract gcd of bases in num_rat
# 2**(1/3)*6**(1/4) -> 2**(1/3+1/4)*3**(1/4)
pnew = {}
i = 0 # steps through num_rat which may grow
while i < len(num_rat):
bi, ei = num_rat[i]
grow = []
for j in range(i + 1, len(num_rat)):
bj, ej = num_rat[j]
g = igcd(bi, bj)
if g != 1:
# 4**r1*6**r2 -> 2**(r1+r2)  *  2**r1 *  3**r2
# this might have a gcd with something else
e = ei + ej
if e.q == 1:
coeff *= Pow(g, e)
else:
if e.p > e.q:
e_i, ep = divmod(e.p, e.q) # change e in place
coeff *= Pow(g, e_i)
e = Rational(ep, e.q)
grow.append((g, e))
# update the jth item
num_rat[j] = (bj//g, ej)
# update bi that we are checking with
bi = bi//g
if bi is S.One:
break
if bi is not S.One:
obj = Pow(bi, ei)
if obj.is_Number:
coeff *= obj
else:
if obj.is_Mul: # 12**(1/2) -> 2*sqrt(3)
c, obj = obj.args # expecting only 2 args
coeff *= c
assert obj.is_Pow
bi, ei = obj.args
pnew.setdefault(ei, []).append(bi)

num_rat.extend(grow)
i += 1

# combine bases of the new powers
for e, b in pnew.iteritems():
pnew[e] = Mul(*b)

# see if there is a base with matching coefficient
# that the -1 can be joined with
if neg1e:
p = Pow(S.NegativeOne, neg1e)
if p.is_Number:
coeff *= p
else:
if p.is_Mul:
c, p = p.args
coeff *= c
assert p.is_Pow and p.base is S.NegativeOne
neg1e = p.args[1]
for e, b in pnew.iteritems():
if e == neg1e and b.is_positive:
pnew[e] = -b
break
else:
c_part.append(p)

# add all the pnew powers
c_part.extend([Pow(b, e) for e, b in pnew.iteritems()])

# oo, -oo
if (coeff is S.Infinity) or (coeff is S.NegativeInfinity):
new_c_part = []
coeff_sign = 1
for t in c_part:
if t.is_positive:
continue
if t.is_negative:
coeff_sign *= -1
continue
new_c_part.append(t)
c_part = new_c_part
new_nc_part = []
for t in nc_part:
if t.is_positive:
continue
if t.is_negative:
coeff_sign *= -1
continue
new_nc_part.append(t)
nc_part = new_nc_part
coeff *= coeff_sign

# zoo
if coeff is S.ComplexInfinity:
# zoo might be
#   unbounded_real + bounded_im
#   bounded_real + unbounded_im
#   unbounded_real + unbounded_im
# and non-zero real or imaginary will not change that status.
c_part = [c for c in c_part if not (c.is_nonzero and
c.is_real is not None)]
nc_part = [c for c in nc_part if not (c.is_nonzero and
c.is_real is not None)]

# 0
elif coeff is S.Zero:
# we know for sure the result will be 0
return [coeff], [], order_symbols

# order commutative part canonically
c_part.sort(key=cmp_to_key(Basic.compare))

# current code expects coeff to be always in slot-0
if coeff is not S.One:
c_part.insert(0, coeff)

# we are done
if len(c_part)==2 and c_part[0].is_Number and c_part[1].is_Add:
# 2*(1+a) -> 2 + 2 * a
coeff = c_part[0]
c_part = [Add(*[coeff*f for f in c_part[1].args])]

return c_part, nc_part, order_symbols

def _eval_power(b, e):
if e.is_Number:
if b.is_commutative:
if e.is_Integer:
# (a*b)**2 -> a**2 * b**2
return Mul(*[s**e for s in b.args])

if e.is_rational:
coeff, rest = b.as_coeff_mul()
rest = list(rest)
unk=[]
nonneg=[]
neg=[]
for bi in rest:
if not bi.is_negative is None: #then we know the sign
if bi.is_negative:
neg.append(bi)
else:
nonneg.append(bi)
else:
unk.append(bi)
if len(unk) == len(rest) or len(neg) == len(rest) == 1:
# if all terms were unknown there is nothing to pull
# out except maybe the coeff; if there is a single
# negative term, this is the base case which cannot
# be processed further
if coeff.is_negative:
coeff *= -1
rest[0] = -rest[0]
if coeff is S.One:
return None
return Mul(Pow(coeff, e), Pow(Mul(*rest), e))

# otherwise return the new expression expanding out the
# known terms; those that are not known can be expanded
# out with separate() but this will introduce a lot of
# "garbage" that is needed to keep one on the same branch
# as the unexpanded expression. The negatives are brought
# out with a negative sign added and a negative left behind
# in the unexpanded terms if there were an odd number of
# negatives.
if neg:
neg = [-w for w in neg]
if coeff.is_negative:
coeff = -coeff
unk.append(S.NegativeOne)
if len(neg) % 2:
unk.append(S.NegativeOne)
return Mul(*[Pow(s, e) for s in nonneg + neg + [coeff]])* \
Pow(Mul(*unk), e)

coeff, rest = b.as_coeff_mul()
if coeff is not S.One:
# (2*a)**3 -> 2**3 * a**3
return Mul(Pow(coeff, e), Mul(*[s**e for s in rest]))
elif e.is_Integer:
coeff, rest = b.as_coeff_mul()
if coeff == S.One:
return # the test below for even exponent needs coeff != 1
else:
return Mul(Pow(coeff, e), Pow(Mul(*rest), e))

c, t = b.as_coeff_mul()
if e.is_even and c.is_Number and c < 0:
return Pow((Mul(-c, Mul(*t))), e)

#if e.has(Wild):
#    return Mul(*[t**e for t in b])

@classmethod
def class_key(cls):
return 3, 0, cls.__name__

def _eval_evalf(self, prec):
return AssocOp._eval_evalf(self, prec).expand()

@cacheit
[docs]    def as_two_terms(self):
"""Return head and tail of self.

This is the most efficient way to get the head and tail of an
expression.

- if you want only the head, use self.args[0];
- if you want to process the arguments of the tail then use
self.as_coef_mul() which gives the head and a tuple containing
the arguments of the tail when treated as a Mul.
- if you want the coefficient when self is treated as an Add

>>> from sympy.abc import x, y
>>> (3*x*y).as_two_terms()
(3, x*y)
"""
args = self.args

if len(args) == 1:
return S.One, self
elif len(args) == 2:
return args

else:
return args[0], self._new_rawargs(*args[1:])

@cacheit
def as_coeff_mul(self, *deps):
if deps:
l1 = []
l2 = []
for f in self.args:
if f.has(*deps):
l2.append(f)
else:
l1.append(f)
return self._new_rawargs(*l1), tuple(l2)
coeff, notrat = self.args[0].as_coeff_mul()
if not coeff is S.One:
return coeff, notrat + self.args[1:]
return S.One, self.args

@staticmethod
def _expandsums(sums):
"""
Helper function for _eval_expand_mul.

sums must be a list of instances of Basic.
"""

L = len(sums)
if L == 1:
return sums[0].args
terms = []
left = Mul._expandsums(sums[:L//2])
right = Mul._expandsums(sums[L//2:])

terms = [Mul(a, b) for a in left for b in right]

def _eval_expand_basic(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_basic'):
newterm = term._eval_expand_basic(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_power_exp(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_power_exp'):
newterm = term._eval_expand_power_exp(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_power_base(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_power_base'):
newterm = term._eval_expand_power_base(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_mul(self, deep=True, **hints):
plain, sums, rewrite = [], [], False
for factor in self.args:
if deep:
term = factor.expand(deep=deep, **hints)
if term != factor:
factor = term
rewrite = True

sums.append(factor)
rewrite = True
else:
if factor.is_commutative:
plain.append(factor)
else:
Wrapper = Basic
sums.append(Wrapper(factor))

if not rewrite:
return self
else:
if sums:
terms = Mul._expandsums(sums)
plain = Mul(*plain)
return Add(*[Mul(plain, term) for term in terms])
else:
return Mul(*plain)

def _eval_expand_multinomial(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_multinomial'):
newterm = term._eval_expand_multinomial(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_log(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_log'):
newterm = term._eval_expand_log(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_complex(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_complex'):
newterm = term._eval_expand_complex(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_trig(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_trig'):
newterm = term._eval_expand_trig(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_expand_func(self, deep=True, **hints):
sargs, terms = self.args, []
for term in sargs:
if hasattr(term, '_eval_expand_func'):
newterm = term._eval_expand_func(deep=deep, **hints)
else:
newterm = term
terms.append(newterm)
return self.func(*terms)

def _eval_derivative(self, s):
terms = list(self.args)
factors = []
for i in xrange(len(terms)):
t = terms[i].diff(s)
if t is S.Zero:
continue
factors.append(Mul(*(terms[:i]+[t]+terms[i+1:])))

def _matches_simple(self, expr, repl_dict):
# handle (w*3).matches('x*5') -> {w: x*5/3}
coeff, terms = self.as_coeff_mul()
if len(terms) == 1:
newexpr = self.__class__._combine_inverse(expr, coeff)
return terms[0].matches(newexpr, repl_dict)
return

def matches(self, expr, repl_dict={}, evaluate=False):
expr = sympify(expr)
if self.is_commutative and expr.is_commutative:
return AssocOp._matches_commutative(self, expr, repl_dict, evaluate)
# todo for commutative parts, until then use the default matches method for non-commutative products
return self._matches(expr, repl_dict, evaluate)

def _matches(self, expr, repl_dict={}, evaluate=False):
if evaluate:
return self.subs(repl_dict).matches(expr, repl_dict)

# weed out negative one prefixes
sign = 1
a, b = self.as_two_terms()
if a is S.NegativeOne:
if b.is_Mul:
sign = -sign
else:
# the remainder, b, is not a Mul anymore
return b.matches(-expr, repl_dict, evaluate)
expr = sympify(expr)
if expr.is_Mul and expr.args[0] is S.NegativeOne:
expr = -expr; sign = -sign

if not expr.is_Mul:
# expr can only match if it matches b and a matches +/- 1
if len(self.args) == 2:
# quickly test for equality
if b == expr:
return a.matches(Rational(sign), repl_dict, evaluate)
# do more expensive match
dd = b.matches(expr, repl_dict, evaluate)
if dd == None:
return None
dd = a.matches(Rational(sign), dd, evaluate)
return dd
return None

d = repl_dict.copy()

# weed out identical terms
pp = list(self.args)
ee = list(expr.args)
for p in self.args:
if p in expr.args:
ee.remove(p)
pp.remove(p)

# only one symbol left in pattern -> match the remaining expression
if len(pp) == 1 and isinstance(pp[0], C.Wild):
if len(ee) == 1:
d[pp[0]] = sign * ee[0]
else:
d[pp[0]] = sign * (type(expr)(*ee))
return d

if len(ee) != len(pp):
return None

i = 0
for p, e in zip(pp, ee):
if i == 0 and sign != 1:
try:
e = sign * e
except TypeError:
return None
d = p.matches(e, d, evaluate=not i)
i += 1
if d is None:
return None
return d

@staticmethod
def _combine_inverse(lhs, rhs):
"""
Returns lhs/rhs, but treats arguments like symbols, so things like
oo/oo return 1, instead of a nan.
"""
if lhs == rhs:
return S.One
def check(l, r):
if l.is_Float and r.is_comparable:
return False
if check(lhs, rhs) or check(rhs, lhs):
return S.One
if lhs.is_Mul and rhs.is_Mul:
a = list(lhs.args)
b = [1]
for x in rhs.args:
if x in a:
a.remove(x)
else:
b.append(x)
return Mul(*a)/Mul(*b)
return lhs/rhs

def as_powers_dict(self):
d = {}
for term in self.args:
b, e = term.as_base_exp()
if b not in d:
d[b] = e
else:
d[b] += e
return d

def as_numer_denom(self):
numers, denoms = [],[]
for t in self.args:
n,d = t.as_numer_denom()
numers.append(n)
denoms.append(d)
return Mul(*numers), Mul(*denoms)

[docs]    def as_coeff_Mul(self):
"""Efficiently extract the coefficient of a product. """
coeff, args = self.args[0], self.args[1:]

if coeff.is_Number:
if len(args) == 1:
return coeff, args[0]
else:
return coeff, self._new_rawargs(*args)
else:
return S.One, self

def as_base_exp(self):
e1 = None
bases = []
nc = 0
for m in self.args:
b, e = m.as_base_exp()
if not b.is_commutative:
nc += 1
if e1 is None:
e1 = e
elif e != e1 or nc > 1:
return self, S.One
bases.append(b)
return Mul(*bases), e1

def _eval_is_polynomial(self, syms):
return all(term._eval_is_polynomial(syms) for term in self.args)

def _eval_is_rational_function(self, syms):
return all(term._eval_is_rational_function(syms) for term in self.args)

_eval_is_bounded = lambda self: self._eval_template_is_attr('is_bounded')
_eval_is_commutative = lambda self: self._eval_template_is_attr('is_commutative')
_eval_is_integer = lambda self: self._eval_template_is_attr('is_integer')
_eval_is_comparable = lambda self: self._eval_template_is_attr('is_comparable')

# I*I -> R,  I*I*I -> -I

def _eval_is_real(self):
im_count = 0
is_neither = False
for t in self.args:
if t.is_imaginary:
im_count += 1
continue
t_real = t.is_real
if t_real:
continue
elif t_real is False:
if is_neither:
return None
else:
is_neither = True
else:
return None
if is_neither:
return False

return (im_count % 2 == 0)

def _eval_is_imaginary(self):
im_count = 0
is_neither = False
for t in self.args:
if t.is_imaginary:
im_count += 1
continue
t_real = t.is_real
if t_real:
continue
elif t_real is False:
if is_neither:
return None
else:
is_neither = True
else:
return None
if is_neither:
return False

return (im_count % 2 == 1)

def _eval_is_irrational(self):
for t in self.args:
a = t.is_irrational
if a:
return True
if a is None:
return
return False

def _eval_is_positive(self):
terms = [t for t in self.args if not t.is_positive]
if not terms:
return True
c = terms[0]
if len(terms)==1:
if c.is_nonpositive:
return False
return
r = self._new_rawargs(*terms[1:])
if c.is_negative and r.is_negative:
return True
if r.is_negative and c.is_negative:
return True
# check for nonpositivity, <=0
if c.is_negative and r.is_nonnegative:
return False
if r.is_negative and c.is_nonnegative:
return False
if c.is_nonnegative and r.is_nonpositive:
return False
if r.is_nonnegative and c.is_nonpositive:
return False

def _eval_is_negative(self):
terms = [t for t in self.args if not t.is_positive]
if not terms:
# all terms are either positive -- 2*Symbol('n', positive=T)
#               or     unknown  -- 2*Symbol('x')
if self.is_positive:
return False
else:
return None
c = terms[0]
if len(terms)==1:
return c.is_negative
r = self._new_rawargs(*terms[1:])
# check for nonnegativity, >=0
if c.is_negative and r.is_nonpositive:
return False
if r.is_negative and c.is_nonpositive:
return False
if c.is_nonpositive and r.is_nonpositive:
return False
if c.is_nonnegative and r.is_nonnegative:
return False

def _eval_is_odd(self):
is_integer = self.is_integer

if is_integer:
r = True
for t in self.args:
if t.is_even:
return False
if t.is_odd is None:
r = None
return r

# !integer -> !odd
elif is_integer == False:
return False

def _eval_is_even(self):
is_integer = self.is_integer

if is_integer:
return fuzzy_not(self._eval_is_odd())

elif is_integer == False:
return False

def _eval_subs(self, old, new):

from sympy import sign
from sympy.simplify.simplify import powdenest

if self == old:
return new

def fallback():
"""Return this value when partial subs has failed."""

return self.__class__(*[s._eval_subs(old, new) for s in
self.args])

def breakup(eq):
"""break up powers assuming (not checking) that eq is a Mul:
b**(Rational*e) -> b**e, Rational
commutatives come back as a dictionary {b**e: Rational}
noncommutatives come back as a list [(b**e, Rational)]
"""

(c, nc) = (dict(), list())
for (i, a) in enumerate(Mul.make_args(eq)):
a = powdenest(a)
(b, e) = a.as_base_exp()
if not e is S.One:
(co, _) = e.as_coeff_mul()
b = Pow(b, e/co)
e = co
if a.is_commutative:
if b in c: # handle I and -1 like things where b, e for I is -1, 1/2
c[b] += e
else:
c[b] = e
else:
nc.append([b, e])
return (c, nc)

def rejoin(b, co):
"""
Put rational back with exponent; in general this is not ok, but
since we took it from the exponent for analysis, it's ok to put
it back.
"""

(b, e) = b.as_base_exp()
return Pow(b, e*co)

def ndiv(a, b):
"""if b divides a in an extractive way (like 1/4 divides 1/2
but not vice versa, and 2/5 does not divide 1/3) then return
the integer number of times it divides, else return 0.
"""

if not b.q % a.q or not a.q % b.q:
return int(a/b)
return 0

if not old.is_Mul:
return fallback()

# handle the leading coefficient and use it to decide if anything
# should even be started; we always know where to find the Rational
# so it's a quick test

co_self = self.args[0]
co_old = old.args[0]
if co_old.is_Rational and co_self.is_Rational:
co_xmul = co_self.extract_multiplicatively(co_old)
elif co_old.is_Rational:
co_xmul = None
else:
co_xmul = True

if not co_xmul:
return fallback()

(c, nc) = breakup(self)
(old_c, old_nc) = breakup(old)

# update the coefficients if we had an extraction

if getattr(co_xmul, 'is_Rational', False):
c.pop(co_self)
c[co_xmul] = S.One
old_c.pop(co_old)

# Do some quick tests to see whether we can succeed:
# 1) check for more non-commutative or 2) commutative terms
# 3) ... unmatched non-commutative bases
# 4) ... unmatched commutative terms
# 5) and finally differences in sign
if len(old_nc) > len(nc) or len(old_c) > len(c) or \
set(_[0] for _ in  old_nc).difference(set(_[0] for _ in nc)) or \
set(old_c).difference(set(c)) or \
any(sign(c[b]) != sign(old_c[b]) for b in old_c):
return fallback()

if not old_c:
cdid = None
else:
rat = []
for (b, old_e) in old_c.items():
c_e = c[b]
rat.append(ndiv(c_e, old_e))
if not rat[-1]:
return fallback()
cdid = min(rat)

if not old_nc:
ncdid = None
for i in range(len(nc)):
nc[i] = rejoin(*nc[i])
else:
ncdid = 0  # number of nc replacements we did
take = len(old_nc)  # how much to look at each time
limit = cdid or S.Infinity # max number that we can take
failed = []  # failed terms will need subs if other terms pass
i = 0
while limit and i + take <= len(nc):
hit = False

# the bases must be equivalent in succession, and
# the powers must be extractively compatible on the
# first and last factor but equal inbetween.

rat = []
for j in range(take):
if nc[i + j][0] != old_nc[j][0]:
break
elif j == 0:
rat.append(ndiv(nc[i + j][1], old_nc[j][1]))
elif j == take - 1:
rat.append(ndiv(nc[i + j][1], old_nc[j][1]))
elif nc[i + j][1] != old_nc[j][1]:
break
else:
rat.append(1)
j += 1
else:
ndo = min(rat)
if ndo:
if take == 1:
if cdid:
ndo = min(cdid, ndo)
nc[i] = Pow(new, ndo)*rejoin(nc[i][0],
nc[i][1] - ndo*old_nc[0][1])
else:
ndo = 1

# the left residual

l = rejoin(nc[i][0], nc[i][1] - ndo*
old_nc[0][1])

# eliminate all middle terms

mid = new

# the right residual (which may be the same as the middle if take == 2)

ir = i + take - 1
r = (nc[ir][0], nc[ir][1] - ndo*
old_nc[-1][1])
if r[1]:
if i + take < len(nc):
nc[i:i + take] = [l*mid, r]
else:
r = rejoin(*r)
nc[i:i + take] = [l*mid*r]
else:

# there was nothing left on the right

nc[i:i + take] = [l*mid]

limit -= ndo
ncdid += ndo
hit = True
if not hit:

# do the subs on this failing factor

failed.append(i)
i += 1
else:

if not ncdid:
return fallback()

# although we didn't fail, certain nc terms may have
# failed so we rebuild them after attempting a partial
# subs on them

failed.extend(range(i, len(nc)))
for i in failed:
nc[i] = rejoin(*nc[i]).subs(old, new)

# rebuild the expression

if cdid is None:
do = ncdid
elif ncdid is None:
do = cdid
else:
do = min(ncdid, cdid)

margs = []
for b in c:
if b in old_c:

# calculate the new exponent

e = c[b] - old_c[b]*do
margs.append(rejoin(b, e))
else:
margs.append(rejoin(b.subs(old, new), c[b]))
if cdid and not ncdid:

# in case we are replacing commutative with non-commutative,
# we want the new term to come at the front just like the
# rest of this routine

margs = [Pow(new, cdid)] + margs
return Mul(*margs)*Mul(*nc)

def _eval_nseries(self, x, n, logx):
from sympy import powsimp
terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
return powsimp(Mul(*terms).expand(), combine='exp', deep=True)

return Mul(*[t.as_leading_term(x) for t in self.args])

def _eval_conjugate(self):
return Mul(*[t.conjugate() for t in self.args])

def _sage_(self):
s = 1
for x in self.args:
s *= x._sage_()
return s

from numbers import Rational, igcd
from power import Pow
from sympify import sympify