Ask Your Question
0

factoring out a term - not only a variable - in a symbolic expression

asked 2022-04-25 12:46:26 +0100

Albert_Zweistein gravatar image

updated 2022-04-25 13:24:29 +0100

Max Alekseyev gravatar image

I know that I can use collect to "factor out" a variable or arrange a symbolic expressions in powers of that variable.

But if there's a term, e.g. x+1, as common factor, I don't know how to factor it out. collect doesn't work for this, e.g.:

sage: var('c0 c1 c2 x0 x1 x')
sage: (c0 * (x+1) + c1*(x+1) + c2*c1*x - c0*c1*(x+1)).collect(x+1)
output: -c0*c1*(x + 1) + c1*c2*x + c0*(x + 1) + c1*(x + 1)

What I would like to have here is:

(-c0*c1 + c0 + c1)*(x + 1) + c1*c2*x
edit retag flag offensive close merge delete

3 Answers

Sort by ยป oldest newest most voted
2

answered 2022-04-25 17:17:22 +0100

Emmanuel Charpentier gravatar image

A tad awkwardly :

sage: sum([v[0]*SR("u")^v[1] for v in (c0 * (x+1) + c1*(x+1) + c2*c1*x - c0*c1*
        (x+1)).subs(x+1==SR("u")).coefficients(SR("u"))]).subs(SR("u")==x+1)
c1*c2*x - (c0*c1 - c0 - c1)*(x + 1)

This can be programmed as :

def mycollect(e1, e2):
    """
    Collects e2 coefficients in e1, even when e2 is not a factor of
    some of the terms of e1
    """
    g = maxima_calculus.gensym()._sage_()
    while g in e1.variables(): g = maxima_calculus.gensym()._sage_()
    return sum([v[0]*g^v[1]
                for v in e1.subs(e2==g).coefficients(g)]).subs(g==e2)

Test :

sage: mycollect(c0 * (x+1) + c1*(x+1) + c2*c1*x - c0*c1*(x+1),x+1)
c1*c2*x - (c0*c1 - c0 - c1)*(x + 1)

This (with some possible generalizations) might become a new method for symbolic expressions. Advice ?

HTH,

edit flag offensive delete link more

Comments

I like your solution. It also works e.g. for this

sage: mycollect(c0 * (cos(x)+exp(x)) + c1*(cos(x)+exp(x)) + c2*c1*x - c0*c1*(cos(x)+exp(x)), cos(x)+exp(x))

c1*c2*x - (c0*c1 - c0 - c1)*(cos(x) + e^x)
Albert_Zweistein gravatar imageAlbert_Zweistein ( 2022-04-26 08:45:27 +0100 )edit

Unfortunately this seems not to work when factoring out multiple factors at the same time:

sage: var('x1 y1 y2 y3')
sage: mycollect(x1*y1*y2 + y3*x1*y1, x1*y1)
x1*y1*y2 + x1*y1*y3

I would like to get

(y2 + y3) * x1*y1
here.

I also tried

sage: mycollect(mycollect(x1*y1*y2 + y3*x1*y1, x1), y1)
(x1*y2 + x1*y3)*y1

Has anyone a solution that also works for this ?

Albert_Zweistein gravatar imageAlbert_Zweistein ( 2022-05-01 14:47:29 +0100 )edit
1

answered 2022-04-25 15:57:33 +0100

Max Alekseyev gravatar image

updated 2022-04-25 16:00:56 +0100

First off, notice the result you request is rather arbitrary - why it's (-c0*c1 + c0 + c1)*(x + 1) + c1*c2*x but not (-c0*c1 + c0 + c1 + c1*c2)*(x + 1) - c1*c2?

The latter can be achieved via using a new variable $y=x+1$:

var('c0 c1 c2 x0 x1 x y')
(c0 * (x+1) + c1*(x+1) + c2*c1*x - c0*c1*(x+1)).subs({x:y-1}).collect(y)

However, if your expressions are polynomial it may be worth to work in polynomial rather than symbolic ring, which provides much advanced functionality with respect to handling polynomials.

edit flag offensive delete link more

Comments

Yes, you're right, it would have been clearer if I had written e.g. -c0c1(x + 1) + c1c2x1 + c0(x + 1) + c1(x + 1). But sometimes I get long results from a calculation and want to factor out only some terms that I see in the output as common factors of some products (not those that could be used when writing the expression in another way).

Albert_Zweistein gravatar imageAlbert_Zweistein ( 2022-04-26 08:39:34 +0100 )edit

I have found a solution that seems to work for me, but the code may be complicated or even buggy. See my updated post.

Albert_Zweistein gravatar imageAlbert_Zweistein ( 2022-10-07 06:30:06 +0100 )edit
0

answered 2022-05-01 16:49:33 +0100

Albert_Zweistein gravatar image

updated 2022-10-08 07:18:14 +0100

I reformulate my wish for a function - call it my_collect2(expr, fact) - more precisely:

Assume expr is a sum of products. Each factor in each product can be a constant, a single variable or any other expression like a+b, xy, sin(axy), ... If the sum contains also a single variable (rather than a product), I consider this variable also as "product", e.g. in (a+b)x + x the last x shall also be considered as "product" and so this whole expression is a sum of products. For this expression

my_collect2((a+b)*x + x, x)

shall yield

(a+b+1) * x

Generally assume that fact is also any expression. If fact occures as common factor in two or more of the products, fact should be factored out by the call my_collect2(expr, fact).

Now assume expr is not a sum of products, but when seeing expr as a tree, there may be one or more sums of products inside. Then the factoring-out shall also be applied recursively to each 'inner' sum of products.

More generally, the recursion shall also be made for a python list of expressions, i.e. each expression in the list shall be searched for possible factoring-out's.

Apart from factoring out I do not want to apply any mathemetical rearrangements of the expression (that keep the expression value unchanged). With this I mean e.g. if you have

var('a b c x y z')
my_collect2( (a+b)*(x+y-1) + a+b + 2*c*(x+y) , x*y)

I do not want to get

(a+b+2*c) * (x+y)

since this would first require to rearrange the expression to (a+b)(x+y) + 2c*(x+y) and then factor out (x+y).

Further more, my_collect2 should also be callable with three arguments:

my_collect2(expr, fact, subs_fact)

should work like my_collect2(expr, fact) except that the each factored out term fact should be substituted by subs_fact, e.g.

sage: var('x1 y1 y2 y3 z1 z2 x1y1')
sage: mycollect2(x1*y1*y2 + y3*x1*y1, x1*y1, x1y1)

Expected: (y2 + y3) * x1y1

It should be possible to write such a function my_collect2(), but since I'm new to sage it will probably be easier to do for someone with more experience. My idea is to recursively scan the expression tree for sums of products and then look if fact is a common factor of at least two of the products.

I give some examples:

sage: var('x1 y1 y2 y3 z1 z2')
sage: my_collect2(x1*y1*y2 + y3*x1*y1, x1*y1)

Expected: (y2 + y3) * x1*y1

sage: var('x1 y1 y2 y3 z1')
sage: my_collect2(sin(exp( (x1+cos(y1))*z1*y2 + y3*(x1+cos(y1))*z1  )), (x1+cos(y1))*z1)
Expected: sin(exp( (y2 + y3) * (x1+cos(y1))*z1  ))

sage: my_collect2([z1*(x1+y1) + z2*(x1+y1), (z1-1)*(x1+y1) + x1 + y1 + z2*(x1+y1)], x1+y1)
Expected: [(z1 + z2)*(x1+y1), (z1-1)*(x1+y1) + x1 + y1 + z2*(x1+y1)]

Edit 10/04/2022:

Some months ago I wrote a code that fulfills partly what I want. But as I'm beginner with sagemath, I'm sure that a lot of things in the code are too complicated, not optimal or even buggy. And the code was "in work", but then I didn't have the time to complete it in a way that I thought it was "finished" in some way. But I think just discarding the code wouldn't be good. So here it is:

from sage.symbolic.operators import add_vararg, mul_vararg

def is_exact_positive_integer_multiple_of(x, test_multiple):
    # returns None if x or test_multiple are not exact (containing floating point numbers).
    # Otherwise, returns test_multiple / x if sagemath can prove (by simplification, if necessary) that there
    # is a positive integer n so that test_multiple == n * x. If there is no such integer or there
    # is one but sagemath cannot prove that an integer n > 0 with test_multiple == n * x exists,
    # is_exact_positive_integer_multiple_of(x, test_multiple) returns None.
    # E.g. is_exact_positive_integer_multiple_of(3, 12*(cos(3/7*pi) - cos(2/7*pi) + cos(1/7*pi)))
    # could return 2 since 12*(cos(3/7*pi) - cos(2/7*pi) + cos(1/7*pi)) == 2*3, but
    # it could also return None since sagemath may not be able to prove that this holds.
    if not SR(x).is_exact() or not SR(test_multiple).is_exact():
        return None
    x = SR(x).simplify_full()
    test_multiple = SR(test_multiple).simplify_full()
    if 0 == x:
        return None
    else:
        if ((test_multiple / x).is_integer() or \
                (SR(x).is_constant() and SR(test_multiple).is_constant() and QQbar(test_multiple / x).is_integer())) \
                and test_multiple >= x:
            return test_multiple / x
        else:
            return None

def is_ratio_gerater_equal_one(x, test):
    # returns None if x or test are not exact (containing floating point numbers).
    # Otherwise, returns test / x if sagemath can prove (by simplification, if necessary)
    # that test / x >= 1. If test / x < 1 or sagemath cannot prove that test / x >= 1,
    # is_ratio_gerater_equal_one(x, test) returns None.
    if not SR(x).is_exact() or not SR(test).is_exact():
        return None
    x = SR(x).simplify_full()
    test = SR(test).simplify_full()
    if x.is_zero():
        return None
    else:
        if x.is_integer() and test.is_integer():
            if (x > 0 and test > x-1) or (x < 0 and test < x+1):
                # You may think that this if-interrogation is redundant due to the following
                # interrogation for real x and test.
                # But it seems that sagemath cannot infer that z >= 1 if you declare the following:
                # var("y z")
                # assume(y, "integer")
                # assume(z, "integer", z>0)
                # Without the preceding interrogation
                #     is_ratio_gerater_equal_one(y, y+z-1)
                # would yield False.
                return test / x
            else:
                return None
        elif test / x >= 1 or (SR(x).is_constant() and SR(test).is_constant() and QQbar(test / x) >= 1):
            return test / x
        else:
            return None

def is_proven_integer(x):
    if not SR(x).is_exact():
        return False
    x = SR(x).simplify_full()
    return bool(x.is_integer() or (SR(x).is_constant() and QQbar(x).is_integer()))

def is_proven_nonnegative(x):
    if not SR(x).is_exact():
        return False
    x = SR(x).simplify_full()
    if is_proven_real(x):
        rr = True
    elif x.operator() is operator.pow and is_proven_real(x.operands()[0]) and is_proven_real(x.operands()[1]) and \
            (is_proven_nonnegative(x.operands()[0]) or is_proven_integer(x.operands()[1])):
        rr = True
    else:
        rr = False
    return rr and bool(x >= 0) or (x.is_integer() and bool(x > -1)) or \
            (SR(x).is_constant() and QQbar(x).is_real() and bool(QQbar(x) >= 0))

def is_proven_real(x):
    if not SR(x).is_exact():
        return False
    x = SR(x).simplify_full()
    return bool(x.is_real() or (SR(x).is_constant() and QQbar(x).is_real())) or \
                (x.operator() is operator.pow and x.operands()[0].is_real() and x.operands()[1].is_real() and \
                (is_proven_nonnegative(x.operands()[0]) or is_proven_integer(x.operands()[1])))

def is_proven_positive(x):
    if not SR(x).is_exact():
        return False
    x = SR(x).simplify_full()
    return bool((is_proven_real(x) and x > 0) or (SR(x).is_constant() and QQbar(x) > 0))

def can_contract_consecutive_exponentiation(x, y, z):
    # returns True if a sufficient condition for (x ** y) ** z == x ** (y*z) could be found, otherwise False.
    return is_proven_integer(z) or (is_proven_nonnegative(x) and is_proven_real(y) and is_proven_real(z))

def try_concatenate_consecutive_exponentiations(expr):
    if expr.operator() is not operator.pow:
        return expr, 1
    else:
        base = expr.operands()[0]
        expn = expr.operands()[1]
        while base.operator() is operator.pow and \
                can_contract_consecutive_exponentiation(base.operands()[0], base.operands()[1], expn):
            expn = base.operands()[1] * expn
            base = base.operands()[0]
        return base, expn

def can_expand_power_of_product(x, y, z):
    # returns True if a sufficient condition for (x*y) ** z == x**z * y**z could be found, otherwise False.
    return is_proven_integer(z) or (is_proven_nonnegative(x) and is_proven_nonnegative(y) and is_proven_real(z))

def remove_power_of_factor_from_product_or_power(expr, factr, exponent, max_exponent = False,
                                                 only_from_integer_powers = False):
    if max_exponent:
        return remove_power_of_factor_from_product_or_power_max_exponent(expr, factr, exponent, only_from_integer_powers)
    else:
        r, ee = remove_power_of_factor_from_product_or_power_no_max_exponent(expr, factr, exponent,
                                                                             only_from_integer_powers)
        return r, None if r is None else 1

def remove_power_of_factor_from_product_or_power_no_max_exponent(expr, factr, exponent, only_from_integer_powers):
    # factr must contain at least one variable, and this variable must also be in expr.
    # If expr == factr ** exponent, this function returns the pair 1, 1 .
    # Otherwise: 'expr' must be a product or an exponentiation (power function).
    # If 'factr' is a product or an exponentiation, look at all factors inside 'factr'
    # that contain a variable (no constant factors), otherwise look at 'factr' as a whole. Under some conditions,
    # remove_power_of_factor_from_product_or_power_no_max_exponent() returns the two values
    #    expr / (factr ** exponent)   ,    nn
    # otherwise, it returns None, None. These conditions are:
    # --- if only_from_integer_powers==True:
    #     all factr ???????????????????????????????????????????????????????????????????

    if 0 == exponent or SR(factr).is_constant() or SR(expr).is_constant():   # SR(factr).is_constant(): if factr is constant, i.e. if it contains no variables
        return None, None

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

    if bool(facop == mul_vararg):
        p = expr
        intpow = None
        for fac in factr.operands():
            if not can_expand_power_of_product(SR(factr / fac).simplify(), fac, exponent):
                return None, None
            if SR(fac).is_constant():  # if fac is constant, i.e. if it contains no variables
                p /= fac ** exponent
            else:
                p, nn = remove_power_of_factor_from_product_or_power_no_max_exponent(p, fac, exponent,
                                                                                     only_from_integer_powers)
                if p is None:
                    return None, None
                if intpow is None:
                    intpow = nn
                elif only_from_integer_powers and nn != intpow:
                    return None, None
        return p,  nn if only_from_integer_powers else 0

    op = expr.operator()
    if op is None:
        if bool(expr == factr ** exponent):
            return 1, 1 if only_from_integer_powers else 0
        else:
            return None, None
    elif op is operator.pow and expr.operands()[0]==factr and \
            (is_exact_positive_integer_multiple_of(exponent, expr.operands()[1]) if only_from_integer_powers \
            else is_ratio_gerater_equal_one(exponent, expr.operands()[1])) is not None:
        # The second return value expr.operands()[1] / exponent is real (not complex, nonreal)
        # because of the tested is_exact_positive_integer_multiple_of() or is_ratio_gerater_equal_one()
        # condition.
        return factr ** (expr.operands()[1] - exponent), \
               expr.operands()[1] / exponent if only_from_integer_powers else 0
    elif bool(op == mul_vararg):
        without = None
        for i,fac in enumerate(expr.operands()):
            if bool(fac == factr ** exponent):
                without = 1
                for j,fac2 in enumerate(expr.operands()):
                    if j != i:
                        without *= fac2
                return without,  1 if only_from_integer_powers else 0
            elif fac.operator() is operator.pow:
                # If x and y are variables, x ** (-y) is converted to (x ** y) ** (-1) in sagemath. So check for consecutive exponentiations.
                base, expn = try_concatenate_consecutive_exponentiations(fac)
                if bool(base == factr) and \
                        (is_exact_positive_integer_multiple_of(exponent, expn) if only_from_integer_powers \
                        else is_ratio_gerater_equal_one(exponent, expn)) is not None:
                    without = factr ** (expn - exponent)
                    for j,fac2 in enumerate(expr.operands()):
                        if j != i:
                            without *= fac2
                    # The second return value fac.operands()[1] / exponent is real (not complex, nonreal)
                    # because of the tested is_exact_positive_integer_multiple_of() or is_ratio_gerater_equal_one()
                    # condition.
                    return without,  fac.operands()[1] / exponent if only_from_integer_powers else 0
        return None, None
    else:
        return None, None

def remove_power_of_factor_from_product_or_power_max_exponent(expr, factr, exponent, only_from_integer_powers):

    if SR(exponent).is_zero() or SR(factr).is_constant() or SR(expr).is_constant():   # SR(factr).is_constant(): if factr is constant, i.e. if it contains no variables
        return None, None

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

    if bool(facop == mul_vararg):
        min_abs_exp = None
        p = expr
        for fac in factr.operands():
            if not can_expand_power_of_product(SR(factr / fac).simplify(), fac, exponent):
                return None, None
            if SR(fac).is_constant():  # if fac is constant, i.e. if it contains no variables
                pass  # needn't divide by constant since p will be discarded and re-calculated below
            else:
                p, nn = remove_power_of_factor_from_product_or_power_max_exponent(p, fac, exponent,
                                                                                  only_from_integer_powers)
                if p is None:
                    return None, None
                if min_abs_exp is None or abs(nn) < min_abs_exp:
                    min_abs_exp = abs(nn)
        p = expr
        for fac in factr.operands():
            if SR(fac).is_constant():  # if fac is constant, i.e. if it contains no variables
                p /= fac ** (min_abs_exp * exponent)
            else:
                p, nn = remove_power_of_factor_from_product_or_power_no_max_exponent(p, fac,
                                                                                     min_abs_exp * exponent,
                                                                                     only_from_integer_powers)
                if p is None:
                    return None, None
        return p, min_abs_exp
    op = expr.operator()
    if op is None:
        if bool(expr == factr ** exponent):
            return 1, 1
        else:
            return None, None
    elif op is operator.pow and expr.operands()[0]==factr and only_from_integer_powers and \
            is_exact_positive_integer_multiple_of(exponent, expr.operands()[1]) is not None:
        # The second return ee is real (not complex, nonreal) because of the tested
        # is_exact_positive_integer_multiple_of() condition.
        ee = is_exact_positive_integer_multiple_of(exponent, expr.operands()[1])
        return 1, ee
    elif op is operator.pow and expr.operands()[0]==factr and not only_from_integer_powers and \
            is_ratio_gerater_equal_one(exponent, expr.operands()[1]) is not None:
        ee = is_ratio_gerater_equal_one(exponent, expr.operands()[1])
        # The second return ee is real (not complex, nonreal) because of the tested
        # is_ratio_gerater_equal_one() condition.
        return factr ** (expr.operands()[1] - exponent * ee), ee
    elif bool(op == mul_vararg):
        without = None
        for i,fac in enumerate(expr.operands()):
            if bool(fac == factr ** exponent):
                without = 1
                for j,fac2 in enumerate(expr.operands()):
                    if j != i:
                        without *= fac2
                return without, 1
            elif fac.operator() is operator.pow:
                # If x and y are variables, x ** (-y) is converted to (x ** y) ** (-1) in sagemath. So check for consecutive exponentiations.
                base, expn = try_concatenate_consecutive_exponentiations(fac)
                if bool(base == factr):
                    if only_from_integer_powers and is_exact_positive_integer_multiple_of(exponent, expn) is not None:
                        ee = is_exact_positive_integer_multiple_of(exponent, expn)
                        without = 1
                        for j,fac2 in enumerate(expr.operands()):
                            if j != i:
                                without *= fac2
                        return without, ee
                    elif not only_from_integer_powers and is_ratio_gerater_equal_one(exponent, expn) is not None:
                        ee = is_ratio_gerater_equal_one(exponent, expn)
                        without = factr ** (expn - exponent * ee)
                        for j,fac2 in enumerate(expr.operands()):
                            if j != i:
                                without *= fac2
                        # The second return ee is real (not complex, nonreal)
                        # because of the tested is_ratio_gerater_equal_one() condition.
                        return without, ee
        return None, None
    else:
        return None, None


def factor_out(expr, factr, *, recursive = True, all_integer_powers = False, only_from_integer_powers = False,
                                replace_factor_by = None):

    if type(expr) == list:
        return [factor_out(ex, factr, recursive = recursive, all_integer_powers = all_integer_powers,
                           only_from_integer_powers = only_from_integer_powers, replace_factor_by = replace_factor_by)
                for ex in expr]
    if type(expr) == tuple:
        return tuple(factor_out(ex, factr, recursive = recursive, all_integer_powers = all_integer_powers,
                                only_from_integer_powers = only_from_integer_powers, replace_factor_by = replace_factor_by)
                     for ex in expr)

    if replace_factor_by is None:
        replace_factor_by = factr

    op = expr.operator()
    if op is None or op != add_vararg:
        if op is None or not recursive:
            return expr
        else:
            return op(*[factor_out(ex, factr, recursive = recursive, all_integer_powers = all_integer_powers,
                                   only_from_integer_powers = only_from_integer_powers,
                                   replace_factor_by = replace_factor_by)
                        for ex in expr.operands()])

    summands = expr.operands()
    factors_removed = []

    if all_integer_powers:
        not_containing = summands
        while True:
            fac_removed = []
            summands = not_containing
            max_exponent = None
            for p in summands:
                op = p.operator()
                if op == mul_vararg or (op is operator.pow and p.operands()[0]==factr and p.operands()[1] >= 1):
                    without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, 1, True,
                                                                                      only_from_integer_powers)
                    if without_factor is not None:
                        if max_exponent is None or ee > max_exponent:
                            max_exponent = ee
            if max_exponent is None:
                break 
            not_containing = []  # must be after the preceeding if .. break statements
            for p in summands:
                op = p.operator()
                if op == mul_vararg:
                    without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, max_exponent,
                                                                                      False, only_from_integer_powers)
                    if without_factor is None:
                        if not recursive:
                            not_containing.append(p)
                        else:
                            not_containing.append(op(*[factor_out(ex, factr, recursive = recursive,
                                                                  all_integer_powers = True,
                                                                  only_from_integer_powers = only_from_integer_powers,
                                                                  replace_factor_by = replace_factor_by)
                                                       for ex in p.operands()]))
                    else:
                        op = without_factor.operator()
                        if op is None or not recursive:
                            fac_removed.append(without_factor)
                        else:
                            fac_removed.append(
                                    op(*[factor_out(ex, factr, recursive = recursive,
                                                    all_integer_powers = all_integer_powers,
                                                    only_from_integer_powers = only_from_integer_powers,
                                                    replace_factor_by = replace_factor_by)
                                         for ex in without_factor.operands()]))
                elif op is operator.pow:
                    without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, max_exponent,
                                                                                      False, only_from_integer_powers)
                    if without_factor is None:
                        not_containing.append(op(*[factor_out(ex, factr, recursive = recursive,
                                                              all_integer_powers = True,
                                                              only_from_integer_powers = only_from_integer_powers,
                                                              replace_factor_by = replace_factor_by)
                                                   for ex in p.operands()]))
                    else:
                        op = without_factor.operator()
                        if op is None or not recursive:
                            fac_removed.append(without_factor)
                        else:
                            fac_removed.append(op(*[factor_out(ex, factr, recursive = recursive,
                                                               all_integer_powers = all_integer_powers,
                                                               only_from_integer_powers = only_from_integer_powers,
                                                               replace_factor_by = replace_factor_by)
                                                    for ex in without_factor.operands()]))
                else:
                    if bool(p == factr):
                        fac_removed.append(1)
                    else:
                        not_containing.append(factor_out(p, factr, recursive = recursive, all_integer_powers = True,
                                                         only_from_integer_powers = only_from_integer_powers,
                                                         replace_factor_by = replace_factor_by))
            if 0==len(fac_removed):
                break
            else:
                factors_removed.append( (fac_removed, max_exponent) )
        return sum([sum(fm[0]) * (replace_factor_by ** fm[1]) for fm in factors_removed]) + sum(not_containing)

    else:    

        not_containing = []
        for p in summands:
            op = p.operator()
            if op == mul_vararg:
                without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, 1, False,
                                                                                  only_from_integer_powers)
                if without_factor is None:
                    if not recursive:
                        not_containing.append(p)
                    else:
                        not_containing.append(op(*[factor_out(ex, factr, recursive = recursive,
                                                              all_integer_powers = False,
                                                              only_from_integer_powers = only_from_integer_powers,
                                                              replace_factor_by = replace_factor_by)
                                                   for ex in p.operands()]))
                else:
                    op = without_factor.operator()
                    if op is None or not recursive:
                        factors_removed.append(without_factor)
                    else:
                        factors_removed.append(op(*[factor_out(ex, factr, recursive = recursive,
                                                               all_integer_powers = False,
                                                               only_from_integer_powers = only_from_integer_powers,
                                                               replace_factor_by = replace_factor_by)
                                                    for ex in without_factor.operands()]))
            elif op is operator.pow:
                without_factor, ee = remove_power_of_factor_from_product_or_power(p, factr, 1, all_integer_powers,
                                                                                  only_from_integer_powers)
                if without_factor is None:
                    not_containing.append(op(*[factor_out(ex, factr, recursive = recursive, all_integer_powers = False,
                                                          only_from_integer_powers = only_from_integer_powers,
                                                          replace_factor_by = replace_factor_by)
                                               for ex in p.operands()]))
                else:
                    op = without_factor.operator()
                    if op is None or not recursive:
                        factors_removed.append(without_factor)
                    else:
                        factors_removed.append(op(*[factor_out(ex, factr, recursive = recursive,
                                                               all_integer_powers = all_integer_powers,
                                                               only_from_integer_powers = only_from_integer_powers,
                                                               replace_factor_by = replace_factor_by)
                                                    for ex in without_factor.operands()]))
            else:
                if bool(p == factr):
                    factors_removed.append(1)
                else:
                    not_containing.append(factor_out(p, factr, recursive = recursive,
                                                     all_integer_powers = all_integer_powers,
                                                     only_from_integer_powers = only_from_integer_powers,
                                                     replace_factor_by = replace_factor_by))

        return sum(factors_removed) * replace_factor_by + sum(not_containing)

I give some examples how to use it:

var('a b c x')
expr = 7*a*x^(3/2) + 3*x^(5/2) + b*x^(5/2) + x^(3/2) - 5*x^(-9/2) - c*x^(-9/2) + a*x^5 + 2*x^5 + 7*x^3 + c*x^3 + c*x^(-3) + 9*x^(-3)
factr = x
print("factor_out(", expr, ", ", factr, ")")
for allint in [False, True]:
    for onlyfromint in [False, True]:
        print("allintpow = ", allint, " onlyfromint = ", onlyfromint, " result = ", factor_out(expr, factr, all_integer_powers=allint, only_from_integer_powers=onlyfromint))

The output is:

factor_out( a*x^5 + 2*x^5 + c*x^3 + b*x^(5/2) + 7*x^3 + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2) ,  x )
allintpow =  False  onlyfromint =  False  result =  (a*x^4 + 2*x^4 + c*x^2 + b*x^(3/2) + 7*x^2 + 7*a*sqrt(x) + 3*x^(3/2) + sqrt(x))*x + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
allintpow =  False  onlyfromint =  True  result =  b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + (a*x^4 + 2*x^4 + c*x^2 + 7*x^2)*x + x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  False  result =  (a + 2)*x^5 + (c + 7)*x^3 + (b + 3)*x^(5/2) + (7*a + 1)*x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)
allintpow =  True  onlyfromint =  True  result =  (a + 2)*x^5 + (c + 7)*x^3 + b*x^(5/2) + 7*a*x^(3/2) + 3*x^(5/2) + x^(3/2) + c/x^3 + 9/x^3 - c/x^(9/2) - 5/x^(9/2)

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

Output: x1*y1*(y2 + y3)

var('a b x')
factor_out((a+b)*x + x, x)
Output:  (a + b + 1)*x

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

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2022-04-25 12:46:26 +0100

Seen: 909 times

Last updated: Oct 08 '22