I reformulate my wish for a function - call it my_collect2(expr, fact) - more precisely:
Assume expr is a sum of products. Each factor in each product can be a constant, a single variable or any other expression like a+b, xy, sin(axy), ... If the sum contains also a single variable (rather than a product), I consider this variable also as "product", e.g. in (a+b)x + x the last x shall also be considered as "product" and so this whole expression is a sum of products. For this expression
my_collect2((a+b)*x + x, x)
shall yield
(a+b+1) * x
Generally assume that fact is also any expression. If fact occures as common factor in two or more of the products, fact should be factored out by the call my_collect2(expr, fact).
Now assume expr is not a sum of products, but when seeing expr as a tree, there may be one or more sums of products inside. Then the factoring-out shall also be applied recursively to each 'inner' sum of products.
More generally, the recursion shall also be made for a python list of expressions, i.e. each expression in the list shall be searched for possible factoring-out's.
Apart from factoring out I do not want to apply any mathemetical rearrangements of the expression (that keep the expression value unchanged). With this I mean e.g. if you have
var('a b c x y z')
my_collect2( (a+b)*(x+y-1) + a+b + 2*c*(x+y) , x*y)
I do not want to get
(a+b+2*c) * (x+y)
since this would first require to rearrange the expression to (a+b)(x+y) + 2c*(x+y) and then factor out (x+y).
Further more, my_collect2 should also be callable with three arguments:
my_collect2(expr, fact, subs_fact)
should work like my_collect2(expr, fact) except that the each factored out term fact should be substituted by subs_fact, e.g.
sage: var('x1 y1 y2 y3 z1 z2 x1y1')
sage: mycollect2(x1*y1*y2 + y3*x1*y1, x1*y1, x1y1)
Expected: (y2 + y3) * x1y1
It should be possible to write such a function my_collect2(), but since I'm new to sage it will probably be easier to do for someone with more experience.
My idea is to recursively scan the expression tree for sums of products and then look if fact is a common factor of at least two of the products.
I give some examples:
sage: var('x1 y1 y2 y3 z1 z2')
sage: my_collect2(x1*y1*y2 + y3*x1*y1, x1*y1)
Expected: (y2 + y3) * x1*y1
sage: var('x1 y1 y2 y3 z1')
sage: my_collect2(sin(exp( (x1+cos(y1))*z1*y2 + y3*(x1+cos(y1))*z1 )), (x1+cos(y1))*z1)
Expected: sin(exp( (y2 + y3) * (x1+cos(y1))*z1 ))
sage: my_collect2([z1*(x1+y1) + z2*(x1+y1), (z1-1)*(x1+y1) + x1 + y1 + z2*(x1+y1)], x1+y1)
Expected: [(z1 + z2)*(x1+y1), (z1-1)*(x1+y1) + x1 + y1 + z2*(x1+y1)]
Edit 10/04/2022:
Some months ago I wrote a code that fulfills partly what I want. But as I'm beginner with sagemath, I'm sure that a lot of things in the code are too complicated, not optimal or even buggy. And the code was "in work", but then I didn't have the time to complete it in a way that I thought it was "finished" in some way. But I think just discarding the code wouldn't be good. So here it is:
from sage.symbolic.operators import add_vararg, mul_vararg
def is_exact_positive_integer_multiple_of(x, test_multiple):
# returns None if x or test_multiple are not exact (containing floating point numbers).
# Otherwise, returns test_multiple / x if sagemath can prove (by simplification, if necessary) that there
# is a positive integer n so that test_multiple == n * x. If there is no such integer or there
# is one but sagemath cannot prove that an integer n > 0 with test_multiple == n * x exists,
# is_exact_positive_integer_multiple_of(x, test_multiple) returns None.
# E.g. is_exact_positive_integer_multiple_of(3, 12*(cos(3/7*pi) - cos(2/7*pi) + cos(1/7*pi)))
# could return 2 since 12*(cos(3/7*pi) - cos(2/7*pi) + cos(1/7*pi)) == 2*3, but
# it could also return None since sagemath may not be able to prove that this holds.
if not SR(x).is_exact() or not SR(test_multiple).is_exact():
return None
x = SR(x).simplify_full()
test_multiple = SR(test_multiple).simplify_full()
if 0 == x:
return None
else:
if ((test_multiple / x).is_integer() or \
(SR(x).is_constant() and SR(test_multiple).is_constant() and QQbar(test_multiple / x).is_integer())) \
and test_multiple >= x:
return test_multiple / x
else:
return None
def is_ratio_gerater_equal_one(x, test):
# returns None if x or test are not exact (containing floating point numbers).
# Otherwise, returns test / x if sagemath can prove (by simplification, if necessary)
# that test / x >= 1. If test / x < 1 or sagemath cannot prove that test / x >= 1,
# is_ratio_gerater_equal_one(x, test) returns None.
if not SR(x).is_exact() or not SR(test).is_exact():
return None
x = SR(x).simplify_full()
test = SR(test).simplify_full()
if x.is_zero():
return None
else:
if x.is_integer() and test.is_integer():
if (x > 0 and test > x-1) or (x < 0 and test < x+1):
# You may think that this if-interrogation is redundant due to the following
# interrogation for real x and test.
# But it seems that sagemath cannot infer that z >= 1 if you declare the following:
# var("y z")
# assume(y, "integer")
# assume(z, "integer", z>0)
# Without the preceding interrogation
# is_ratio_gerater_equal_one(y, y+z-1)
# would yield False.
return test / x
else:
return None
elif test / x >= 1 or (SR(x).is_constant() and SR(test).is_constant() and QQbar(test / x) >= 1):
return test / x
else:
return None
def is_proven_integer(x):
if not SR(x).is_exact():
return False
x = SR(x).simplify_full()
return bool(x.is_integer() or (SR(x).is_constant() and QQbar(x).is_integer()))
def is_proven_nonnegative(x):
if not SR(x).is_exact():
return False
x = SR(x).simplify_full()
if is_proven_real(x):
rr = True
elif x.operator() is operator.pow and is_proven_real(x.operands()[0]) and is_proven_real(x.operands()[1]) and \
(is_proven_nonnegative(x.operands()[0]) or is_proven_integer(x.operands()[1])):
rr = True
else:
rr = False
return rr and bool(x >= 0) or (x.is_integer() and bool(x > -1)) or \
(SR(x).is_constant() and QQbar(x).is_real() and bool(QQbar(x) >= 0))
def is_proven_real(x):
if not SR(x).is_exact():
return False
x = SR(x).simplify_full()
return bool(x.is_real() or (SR(x).is_constant() and QQbar(x).is_real())) or \
(x.operator() is operator.pow and x.operands()[0].is_real() and x.operands()[1].is_real() and \
(is_proven_nonnegative(x.operands()[0]) or is_proven_integer(x.operands()[1])))
def is_proven_positive(x):
if not SR(x).is_exact():
return False
x = SR(x).simplify_full()
return bool((is_proven_real(x) and x > 0) or (SR(x).is_constant() and QQbar(x) > 0))
def can_contract_consecutive_exponentiation(x, y, z):
# returns True if a sufficient condition for (x ** y) ** z == x ** (y*z) could be found, otherwise False.
return is_proven_integer(z) or (is_proven_nonnegative(x) and is_proven_real(y) and is_proven_real(z))
def try_concatenate_consecutive_exponentiations(expr):
if expr.operator() is not operator.pow:
return expr, 1
else:
base = expr.operands()[0]
expn = expr.operands()[1]
while base.operator() is operator.pow and \
can_contract_consecutive_exponentiation(base.operands()[0], base.operands()[1], expn):
expn = base.operands()[1] * expn
base = base.operands()[0]
return base, expn
def can_expand_power_of_product(x, y, z):
# returns True if a sufficient condition for (x*y) ** z == x**z * y**z could be found, otherwise False.
return is_proven_integer(z) or (is_proven_nonnegative(x) and is_proven_nonnegative(y) and is_proven_real(z))
def remove_power_of_factor_from_product_or_power(expr, factr, exponent, max_exponent = False,
only_from_integer_powers = False):
if max_exponent:
return remove_power_of_factor_from_product_or_power_max_exponent(expr, factr, exponent, only_from_integer_powers)
else:
r, ee = remove_power_of_factor_from_product_or_power_no_max_exponent(expr, factr, exponent,
only_from_integer_powers)
return r, None if r is None else 1
def remove_power_of_factor_from_product_or_power_no_max_exponent(expr, factr, exponent, only_from_integer_powers):
# factr must contain at least one variable, and this variable must also be in expr.
# If expr == factr ** exponent, this function returns the pair 1, 1 .
# Otherwise: 'expr' must be a product or an exponentiation (power function).
# If 'factr' is a product or an exponentiation, look at all factors inside 'factr'
# that contain a variable (no constant factors), otherwise look at 'factr' as a whole. Under some conditions,
# remove_power_of_factor_from_product_or_power_no_max_exponent() returns the two values
# expr / (factr ** exponent) , nn
# otherwise, it returns None, None. These conditions are:
# --- if only_from_integer_powers==True:
# all factr ???????????????????????????????????????????????????????????????????
if 0 == exponent or SR(factr).is_constant() or SR(expr).is_constant(): # SR(factr).is_constant(): if factr is constant, i.e. if it contains no variables
return None, None
facop = factr.operator()
if facop is operator.pow:
factr, exponent = try_concatenate_consecutive_exponentiations(factr ** exponent)
facop = factr.operator()
if bool(facop == mul_vararg):
p = expr
intpow = None
for fac in factr.operands():
if not can_expand_power_of_product(SR(factr / fac).simplify(), fac, exponent):
return None, None
if SR(fac).is_constant(): # if fac is constant, i.e. if it contains no variables
p /= fac ** exponent
else:
p, nn = remove_power_of_factor_from_product_or_power_no_max_exponent(p, fac, exponent,
only_from_integer_powers)
if p is None:
return None, None
if intpow is None:
intpow = nn
elif only_from_integer_powers and nn != intpow:
return None, None
return p, nn if only_from_integer_powers else 0
op = expr.operator()
if op is None:
if bool(expr == factr ** exponent):
return 1, 1 if only_from_integer_powers else 0
else:
return None, None
elif op is operator.pow and expr.operands()[0]==factr and \
(is_exact_positive_integer_multiple_of(exponent, expr.operands()[1]) if only_from_integer_powers \
else is_ratio_gerater_equal_one(exponent, expr.operands()[1])) is not None:
# The second return value expr.operands()[1] / exponent is real (not complex, nonreal)
# because of the tested is_exact_positive_integer_multiple_of() or is_ratio_gerater_equal_one()
# condition.
return factr ** (expr.operands()[1] - exponent), \
expr.operands()[1] / exponent if only_from_integer_powers else 0
elif bool(op == mul_vararg):
without = None
for i,fac in enumerate(expr.operands()):
if bool(fac == factr ** exponent):
without = 1
for j,fac2 in enumerate(expr.operands()):
if j != i:
without *= fac2
return without, 1 if only_from_integer_powers else 0
elif fac.operator() is operator.pow:
# If x and y are variables, x ** (-y) is converted to (x ** y) ** (-1) in sagemath. So check for consecutive exponentiations.
base, expn = try_concatenate_consecutive_exponentiations(fac)
if bool(base == factr) and \
(is_exact_positive_integer_multiple_of(exponent, expn) if only_from_integer_powers \
else is_ratio_gerater_equal_one(exponent, expn)) is not None:
without = factr ** (expn - exponent)
for j,fac2 in enumerate(expr.operands()):
if j != i:
without *= fac2
# The second return value fac.operands()[1] / exponent is real (not complex, nonreal)
# because of the tested is_exact_positive_integer_multiple_of() or is_ratio_gerater_equal_one()
# condition.
return without, fac.operands()[1] / exponent if only_from_integer_powers else 0
return None, None
else:
return None, None
def remove_power_of_factor_from_product_or_power_max_exponent(expr, factr, exponent, only_from_integer_powers):
if SR(exponent).is_zero() or SR(factr).is_constant() or SR(expr).is_constant(): # SR(factr).is_constant(): if factr is constant, i.e. if it contains no variables
return None, None
facop = factr.operator()
if facop is operator.pow:
factr, exponent = try_concatenate_consecutive_exponentiations(factr ** exponent)
facop = factr.operator()
if bool(facop == mul_vararg):
min_abs_exp = None
p = expr
for fac in factr.operands():
if not can_expand_power_of_product(SR(factr / fac).simplify(), fac, exponent):
return None, None
if SR(fac).is_constant(): # if fac is constant, i.e. if it contains no variables
pass # needn't divide by constant since p will be discarded and re-calculated below
else:
p, nn = remove_power_of_factor_from_product_or_power_max_exponent(p, fac, exponent,
only_from_integer_powers)
if p is None:
return None, None
if min_abs_exp is None or abs(nn) < min_abs_exp:
min_abs_exp = abs(nn)
p = expr
for fac in factr.operands():
if SR(fac).is_constant(): # if fac is constant, i.e. if it contains no variables
p /= fac ** (min_abs_exp * exponent)
else:
p, nn = remove_power_of_factor_from_product_or_power_no_max_exponent(p, fac,
min_abs_exp * exponent,
only_from_integer_powers)
if p is None:
return None, None
return p, min_abs_exp
op = expr.operator()
if op is None:
if bool(expr == factr ** exponent):
return 1, 1
else:
return None, None
elif op is operator.pow and expr.operands()[0]==factr and only_from_integer_powers and \
is_exact_positive_integer_multiple_of(exponent, expr.operands()[1]) is not None:
# The second return ee is real (not complex, nonreal) because of the tested
# is_exact_positive_integer_multiple_of() condition.
ee = is_exact_positive_integer_multiple_of(exponent, expr.operands()[1])
return 1, ee
elif op is operator.pow and expr.operands()[0]==factr and not only_from_integer_powers and \
is_ratio_gerater_equal_one(exponent, expr.operands()[1]) is not None:
ee = is_ratio_gerater_equal_one(exponent, expr.operands()[1])
# The second return ee is real (not complex, nonreal) because of the tested
# is_ratio_gerater_equal_one() condition.
return factr ** (expr.operands()[1] - exponent * ee), ee
elif bool(op == mul_vararg):
without = None
for i,fac in enumerate(expr.operands()):
if bool(fac == factr ** exponent):
without = 1
for j,fac2 in enumerate(expr.operands()):
if j != i:
without *= fac2
return without, 1
elif fac.operator() is operator.pow:
# If x and y are variables, x ** (-y) is converted to (x ** y) ** (-1) in sagemath. So check for consecutive exponentiations.
base, expn = try_concatenate_consecutive_exponentiations(fac)
if bool(base == factr):
if only_from_integer_powers and is_exact_positive_integer_multiple_of(exponent, expn) is not None:
ee = is_exact_positive_integer_multiple_of(exponent, expn)
without = 1
for j,fac2 in enumerate(expr.operands()):
if j != i:
without *= fac2
return without, ee
elif not only_from_integer_powers and is_ratio_gerater_equal_one(exponent, expn) is not None:
ee = is_ratio_gerater_equal_one(exponent, expn)
without = factr ** (expn - exponent * ee)
for j,fac2 in enumerate(expr.operands()):
if j != i:
without *= fac2
# The second return ee is real (not complex, nonreal)
# because of the tested is_ratio_gerater_equal_one() condition.
return without, ee
return None, None
else:
return None, None
def factor_out(expr, factr, *, recursive = True, all_integer_powers = False, only_from_integer_powers = False,
replace_factor_by = None):
if type(expr) == list:
return [factor_out(ex, factr, recursive = recursive, all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers, replace_factor_by = replace_factor_by)
for ex in expr]
if type(expr) == tuple:
return tuple(factor_out(ex, factr, recursive = recursive, all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers, replace_factor_by = replace_factor_by)
for ex in expr)
if replace_factor_by is None:
replace_factor_by = factr
op = expr.operator()
if op is None or op != add_vararg:
if op is None or not recursive:
return expr
else:
return op(*[factor_out(ex, factr, recursive = recursive, all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in expr.operands()])
summands = expr.operands()
factors_removed = []
if all_integer_powers:
not_containing = summands
while True:
fac_removed = []
summands = not_containing
max_exponent = None
for p in summands:
op = p.operator()
if op == mul_vararg or (op is operator.pow and p.operands()[0]==factr and p.operands()[1] >= 1):
without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, 1, True,
only_from_integer_powers)
if without_factor is not None:
if max_exponent is None or ee > max_exponent:
max_exponent = ee
if max_exponent is None:
break
not_containing = [] # must be after the preceeding if .. break statements
for p in summands:
op = p.operator()
if op == mul_vararg:
without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, max_exponent,
False, only_from_integer_powers)
if without_factor is None:
if not recursive:
not_containing.append(p)
else:
not_containing.append(op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = True,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in p.operands()]))
else:
op = without_factor.operator()
if op is None or not recursive:
fac_removed.append(without_factor)
else:
fac_removed.append(
op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in without_factor.operands()]))
elif op is operator.pow:
without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, max_exponent,
False, only_from_integer_powers)
if without_factor is None:
not_containing.append(op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = True,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in p.operands()]))
else:
op = without_factor.operator()
if op is None or not recursive:
fac_removed.append(without_factor)
else:
fac_removed.append(op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in without_factor.operands()]))
else:
if bool(p == factr):
fac_removed.append(1)
else:
not_containing.append(factor_out(p, factr, recursive = recursive, all_integer_powers = True,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by))
if 0==len(fac_removed):
break
else:
factors_removed.append( (fac_removed, max_exponent) )
return sum([sum(fm[0]) * (replace_factor_by ** fm[1]) for fm in factors_removed]) + sum(not_containing)
else:
not_containing = []
for p in summands:
op = p.operator()
if op == mul_vararg:
without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, 1, False,
only_from_integer_powers)
if without_factor is None:
if not recursive:
not_containing.append(p)
else:
not_containing.append(op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = False,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in p.operands()]))
else:
op = without_factor.operator()
if op is None or not recursive:
factors_removed.append(without_factor)
else:
factors_removed.append(op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = False,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in without_factor.operands()]))
elif op is operator.pow:
without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, 1, all_integer_powers,
only_from_integer_powers)
if without_factor is None:
not_containing.append(op(*[factor_out(ex, factr, recursive = recursive, all_integer_powers = False,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in p.operands()]))
else:
op = without_factor.operator()
if op is None or not recursive:
factors_removed.append(without_factor)
else:
factors_removed.append(op(*[factor_out(ex, factr, recursive = recursive,
all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by)
for ex in without_factor.operands()]))
else:
if bool(p == factr):
factors_removed.append(1)
else:
not_containing.append(factor_out(p, factr, recursive = recursive,
all_integer_powers = all_integer_powers,
only_from_integer_powers = only_from_integer_powers,
replace_factor_by = replace_factor_by))
return sum(factors_removed) * replace_factor_by + sum(not_containing)
I give some examples how to use it:
var('a b c x')
expr = 7*a*x^(3/2) + 3*x^(5/2) + b*x^(5/2) + x^(3/2) - 5*x^(-9/2) - c*x^(-9/2) + a*x^5 + 2*x^5 + 7*x^3 + c*x^3 + c*x^(-3) + 9*x^(-3)
factr = x
print("factor_out(", expr, ", ", factr, ")")
for allint in [False, True]:
for onlyfromint in [False, True]:
print("allintpow = ", allint, " onlyfromint = ", onlyfromint, " result = ", factor_out(expr, factr, all_integer_powers=allint, only_from_integer_powers=onlyfromint))
The output is:
factor_out( a*x^5 + 2*x^5 + c*x^3 + b*x^(5/2) + 7*x^3 + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2) , x )
allintpow = False onlyfromint = False result = (a*x^4 + 2*x^4 + c*x^2 + b*x^(3/2) + 7*x^2 + 7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
allintpow = False onlyfromint = True result = b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + (a*x^4 + 2*x^4 + c*x^2 + 7*x^2)*x + x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
allintpow = True onlyfromint = False result = (a + 2)*x^5 + (c + 7)*x^3 + (b + 3)*x^(5/2) + (7*a + 1)*x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
allintpow = True onlyfromint = True result = (a + 2)*x^5 + (c + 7)*x^3 + b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
var('x1 y1 y2 y3 z1 z2')
factor_out(x1*y1*y2 + y3*x1*y1, x1*y1)
Output: x1*y1*(y2 + y3)
var('a b x')
factor_out((a+b)*x + x, x)
Output: (a + b + 1)*x
If you find any bugs or have other suggestions, let me know.
