Revision history [back]

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 0 == x:
         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 ((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 \
             ((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 ((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 ((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, factor, 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, factor, factr, exponent, only_from_integer_powers)
     else:
         r, ee = remove_power_of_factor_from_product_or_power_no_max_exponent(expr, factor, 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, factor, factr, exponent, only_from_integer_powers):
     # factor factr must contain at least one variable, and this variable must also be in expr.
     # If expr == factor factr ** exponent, this function returns the pair 1, 1 .
     # Otherwise: 'expr' must be a product or an exponentiation (power function).
     # If 'factor' 'factr' is a product or an exponentiation, look at all factors inside 'factor'
    'factr'
    # that contain a variable (no constant factors), otherwise look at 'factor' 'factr' as a whole. Under some conditions,
     # remove_power_of_factor_from_product_or_power_no_max_exponent() returns the two values
     #    expr / (factor (factr ** exponent)   ,    nn
     # otherwise, it returns None, None. These conditions are:
     # tbd
--- if only_from_integer_powers==True:
    #     all factr ???????????????????????????????????????????????????????????????????
      if 0 == exponent or SR(factor).is_constant() SR(factr).is_constant() or SR(expr).is_constant():   # SR(factor).is_constant(): SR(factr).is_constant(): if factor factr is constant, i.e. if it contains no variables
         return None, None
     if factor.is_zero(): factr.is_zero():  # this test should normally not be necessary due to above test "factor "factr in CC"
         return None, None

     facop = factor.operator()
    factr.operator()
    if facop is operator.pow:
            factor, factr, exponent = try_concatenate_consecutive_exponentiations(factor try_concatenate_consecutive_exponentiations(factr ** exponent)
         facop = factor.operator()
factr.operator()
      if facop == mul_vararg:
         p = expr
         intpow = None
         for fac in factor.operands():
    factr.operands():
            if not can_expand_power_of_product((factor can_expand_power_of_product((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 expr == factor 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]==factor 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 factor factr ** (expr.operands()[1] - exponent), \
                expr.operands()[1] / exponent if only_from_integer_powers else 0
     elif op == mul_vararg:
         without = None
         for i,fac in enumerate(expr.operands()):
             if fac == factor 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 base == factor 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 = factor 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, factor, factr, exponent, only_from_integer_powers):

     if 0 == exponent or SR(factor).is_constant() SR(factr).is_constant() or SR(expr).is_constant():   # SR(factor).is_constant(): SR(factr).is_constant(): if factor factr is constant, i.e. if it contains no variables
         return None, None
     if factor.is_zero(): factr.is_zero():  # this test should normally not be necessary due to above test "factor "factr in CC"
         return None, None

     facop = factor.operator()
    factr.operator()
    if facop is operator.pow:
            factor, factr, exponent = try_concatenate_consecutive_exponentiations(factor try_concatenate_consecutive_exponentiations(factr ** exponent)
         facop = factor.operator()
factr.operator()
      if facop == mul_vararg:
         min_abs_exp = None
         p = expr
         for fac in factor.operands():
    factr.operands():
            if not can_expand_power_of_product((factor can_expand_power_of_product((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 factor.operands():
    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 expr == factor factr ** exponent:
             return 1, 1
         else:
             return None, None
     elif op is operator.pow and expr.operands()[0]==factor 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]==factor 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 factor factr ** (expr.operands()[1] - exponent * ee), ee
     elif op == mul_vararg:
         without = None
         for i,fac in enumerate(expr.operands()):
             if fac == factor 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 base == factor:
    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 = factor 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, factor, factr, *, recursive = True, all_integer_powers = False, only_from_integer_powers = False,
                                 replace_factor_by = None):

     if type(expr) == list:
         return [factor_out(ex, factor, 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, factor, 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 = factor
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, factor, 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]==factor p.operands()[0]==factr and p.operands()[1] >= 1):
                     without_factor, ee = remove_power_of_factor_from_product_or_power(p, factor, 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, factor, 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, factor, 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, factor, 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, factor, factr, max_exponent,
                                                                                       False, only_from_integer_powers)
                     if without_factor is None:
                         not_containing.append(op(*[factor_out(ex, factor, 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, factor, 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 p == factor:
    factr:
                        fac_removed.append(1)
                     else:
                         not_containing.append(factor_out(p, factor, 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, factor, 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, factor, 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, factor, 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, factor, factr, 1, all_integer_powers,
                                                                                   only_from_integer_powers)
                 if without_factor is None:
                     not_containing.append(op(*[factor_out(ex, factor, 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, factor, 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 p == factor:
    factr:
                    factors_removed.append(1)
                 else:
                     not_containing.append(factor_out(p, factor, 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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

If you find any bugs or have other suggestions, let me know.

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 0 == x:
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 ((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 \
            ((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 ((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 ((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
    if factr.is_zero():  # this test should normally not be necessary due to above test "factr in CC"
        return None, None

    facop = factr.operator()
    if facop is operator.pow:
        factr, exponent = try_concatenate_consecutive_exponentiations(factr ** exponent)
        facop = factr.operator()

    if facop bool(facop == mul_vararg:
mul_vararg):
        p = expr
        intpow = None
        for fac in factr.operands():
            if not can_expand_power_of_product((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 expr bool(expr == factr ** exponent:
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 op bool(op == mul_vararg:
mul_vararg):
        without = None
        for i,fac in enumerate(expr.operands()):
            if fac bool(fac == factr ** exponent:
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 base bool(base == factr 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 0 == exponent 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
    if factr.is_zero():  # this test should normally not be necessary due to above test "factr in CC"
        return None, None

    facop = factr.operator()
    if facop is operator.pow:
        factr, exponent = try_concatenate_consecutive_exponentiations(factr ** exponent)
        facop = factr.operator()

    if facop bool(facop == mul_vararg:
mul_vararg):
        min_abs_exp = None
        p = expr
        for fac in factr.operands():
            if not can_expand_power_of_product((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 expr bool(expr == factr ** exponent:
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 op bool(op == mul_vararg:
mul_vararg):
        without = None
        for i,fac in enumerate(expr.operands()):
            if fac bool(fac == factr ** exponent:
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 base bool(base == factr:
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 p bool(p == factr:
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 p bool(p == factr:
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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

If you find any bugs or have other suggestions, let me know.

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 ((x).is_constant() (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 \
            ((x).is_constant() (SR(x).is_constant() and QQbar(x).is_real() and bool(QQbar(x) 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 ((x).is_constant() (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 ((x).is_constant() (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 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((factr 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 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((factr 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 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
= 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
= 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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

If you find any bugs or have other suggestions, let me know.

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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

var('x1 y1 y2 y3 z1 z2') factor_out(x1y1y2 + y3x1y1, x1*y1)

Output: x1y1(y2 + y3)

If you find any bugs or have other suggestions, let me know.

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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

var('x1 y1 y2 y3 z1 z2') factor_out(x1y1y2 + y3x1y1, x1*y1)

Output: x1y1(y2 + y3)

If you find any bugs or have other suggestions, let me know.

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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

var('x1 y1 y2 y3 z1 z2') factor_out(x1y1y2 + y3x1y1, x1*y1)

Output: x1y1(y2 + y3)

If you find any bugs or have other suggestions, let me know.

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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)

var('x1 y1 y2 y3 z1 z2')
factor_out(x1y1y2 factor_out(x1*y1*y2 + y3x1y1, x1*y1)
y3*x1*y1, x1*y1)
  
Output: x1y1(y2 x1*y1*(y2 + y3)y3)

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 d')
expr = x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) - 5*x^(-9/2)
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( 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (7*a + 1)*x^(3/2) + 3*x^(5/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  7*a*x^(3/2) + 3*x^(5/2) + x^(3/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.

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 d')
x')
expr = x^(3/2) 7*a*x^(3/2) + 3*x^(5/2) + 7*a*x^(3/2) b*x^(5/2) + x^(3/2) - 5*x^(-9/2)
5*x^(-9/2) - c*x^(-9/2)
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( b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - c/x^(9/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (7*a*sqrt(x) (b*x^(3/2) + 7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x - 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) + x^(3/2) - c/x^(9/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (b + 3)*x^(5/2) + (7*a + 1)*x^(3/2) + 3*x^(5/2) - c/x^(9/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) - 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.

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)
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^(9/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (b*x^(3/2) + 7*a*sqrt(x) c/x^3 + 3*x^(3/2) + sqrt(x))*x - 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) + x^(3/2) - c/x^(9/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (b + 3)*x^(5/2) + (7*a + 1)*x^(3/2) - c/x^(9/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) 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.

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)
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
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
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
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 = 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 = 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.