How to speed up the expansion of very large polynomials in sympy?

I am working with very large polynomials in sympy and I need to have them in expanded form to find certain terms and coefficients. However, the expansion of these polynomials takes a long time. Is there a fast way to expand polynomials or get certain terms and coefficients in a different way?

I can find the terms in an expanded polynomial fine, but the time to expand it is the limiting factor.

The polynomials are very large, an example might be:
(x + y + z + a + b + c) ** 24

I've tried both sympy.expand(), and Add.as_poly(). And have found Add.as_poly() to be faster but still very slow.

my_poly = (x + y + z + a + b + c) ** 24
# expand using Add.as_poly()
my_poly.as_poly()
# this takes multiple minutes to execute

I'd like to be able to search through the terms in the expanded polynomial to find ones that contain other terms:
(pseudocode) is x**3*yza**2 contained in 500*x**5*y**2*z*a**4*b*c**2
and if it is contained I want to retrieve the coefficient of that term.

I'm looking to either speed up the expansion, or use a different method to find the desired terms in less time.

The coefficients of a multinomial can be determined without actually expanding the expression. (a + b + c + x + y + z)**5, for example, will have terms whose exponents sum to 5. The example you gave must be a term when the coefficient is 15 (but the coefficient would have been 113513400).

def mcoeff(*args):
    """return multinomial coefficient of a terms in the expansion of
    (x_1 + x_2 + x_3 + ... x_n)**m from the given exponents in the
    term of interest

    EXAMPLES
    ========

    >>> from sympy import expand
    >>> from sympy.abc import x, y, z
    >>> eq = (x + y + z)**4
    >>> eq.expand().coeff(x**2*y**2)
    6
    >>> mcoeff(2,2)
    6
    >>> 
    """
    from sympy import binomial
    from sympy.utilities.misc import as_int
    assert all(as_int(i)>=0 for i in args)
    if len(args) == 2:
        return binomial(sum(args), sorted(args)[0])
    if len(set(args)) == 1:
        return factorial(sum(args))/factorial(args[0])**len(args)
    def runsum(a):
        t = 0
        for i in a:
            t += i
            yield t
    return Mul(*[
        binomial(j, k) for j, k in zip(runsum(args), args)])

If you want to find the coefficient including symbol of something like F = x**3*yza**2...well that's harder but not impossible. If you were seeking those factors in the expansion of (a + b + c + x + y + z)**8 (and I chose 8 so there would be more than one term having that factor; there would be only one such term if the exponent were 7 since 3 + 1 + 1 + 2 = sum of exponents of the terms = 7) then that factor would appear with any other factors of x,y,z,a,b,c that would give a total exponent of 8. In this case it is easy since we need to add only 1 so potential terms having factor F are x*F, y*F, z*F, a*F, b*F, c*F. It will be easier to consider the tuples of the exponents for x,y,z,a,b,c which in F are (3,1,1,2,0,0):

x*F -> (4,1,1,2,0,0)
y*F -> (3,2,1,2,0,0)
z*F -> (3,1,2,2,0,0)
a*F -> (3,1,1,3,0,0)
b*F -> (3,1,1,2,1,0)
c*F -> (3,1,1,2,0,1)

each of those will have a numerical coefficient and a symbolic residual (the factor multiplying F) so the coefficient of F in the expansion will be the sum of those:

>>> from sympy.abc import x, y, z, a, b, c, e
>>> coeff = []; t = [3,1,1,2,0,0]
for i,f in enumerate((x, y, z, a, b, c)):
...     t[i] += 1
...     coeff.append(f*mcoeff(*t))
...     t[i] -= 1
>>> Add(*coeff)
1120*a + 3360*b + 3360*c + 840*x + 1680*y + 1680*z
>>> F = x**3*y*z*a**2
>>> ((x + y + z + a + b + c)**8).expand().subs(F, e).coeff(e)  # wait for it...
1120*a + 3360*b + 3360*c + 840*x + 1680*y + 1680*z

It all gets more complicated when you want to find the terms with smaller factors since then you have to distribute the exponent sum deficiency (in this case only 1 short of 8) over all the symbols, but the idea is the same.

(mcoeff was taken from here with some modification)

Recommended topics

Hot tags