better support for series of generating function

[smichr]

I am not aware of a special method in SymPy to expand 1/(1 - x**e1)/(1-x**e2)/.../(1-x**en)/(1-x**(a*e1*e2*...*en). Series does not detect this special case and times increase dramatically as the number of terms requested is increased. Something like the following might be added to handle this case. This would make working with a generating function in combinatorics much more doable (e.g. in the classroom):

def gseries(*e, n=100):
    """return series expansion to order n at 0 of 1/Mul(*[(1-x**i) for
    i in (e1, e2, e3, ..., G)]) where ei | G.

    Examples
    ========

    >>> from sympy import series, Mul
    >>> from sympy.abc import x
    >>> Mul(*[1 - x**i for i in (1, 2, 6)])
    (1 - x)*(1 - x**2)*(1 - x**6)
    >>> series(1/_, n=8)
    1 - x - x**2 + x**3 - x**6 + x**7 + O(x**8)
    >>> _ - gseries(1, 2, 6, n=8)
    O(x**8)
    """
    from sympy.abc import x
    from sympy import Poly, Mul, Add, quo, expand
    k = max(e)
    assert not any(k%i for i in e)
    d = [1 - x**i for i in sorted(e)]
    lcm = d[-1]
    m = len(e)
    num = Poly(Mul(*[quo(lcm, di) for di in d])).as_expr()
    den = Add(*[binomial(m-1+i,i)*x**(k*i) for i in range(n//k + 1)])
    a = expand(num*den).args
    ex = Add(*[i for i in a if i.is_Number or i.as_independent(x)[1].as_base_exp()[1] < n])
    return ex

Here were some timing results for the making change generator:

gseries(1,5,10,25,50,n=n)
n   time/s
10, 1.29
20, 0.04
30, 0.08
40, 0.07
50, 0.56
60, 0.18
70, 0.1
80, 0.08
90, 0.05

series(1/Mul(*[(1-x**i) for i in (1,5,10,25,50)]),n=n)
10, 0.54
20, 0.74
30, 1.79
40, 4.18
50, 8.33
60, 25.58
70, 36.88

[smichr]

Now I am aware: ring series.

def rpoly(e,c=1,x='x',dom=ZZ):
    """return poly in dom(x) from provided exponents; unless
    specified, all coefficients are 1

    Examples
    ========

    >>> rpoly((0,2,5))
    x**5 + x**2 + 1
    >>> rpoly((0,2,5),(2,3,4))
    4*x**5 + 3*x**2 + 2
    """
    from sympy.polys.rings import ring
    R, x = ring('x', dom)
    if c==1:
        c = [1]*len(e)
    return sum(c[i]*x**e[i] for i in range(len(e)))

def rprod(f, prec=None):
    """return product of ring polynomials

    Examples
    ========

    >>> from sympy.polys.rings import ring
    >>> from sympy import ZZ
    >>> _, x = ring('x', ZZ)
    >>> rprod((1 - x, 1 - x**2, 1 - x**5))
    -x**8 + x**7 + x**6 - x**5 + x**3 - x**2 - x + 1
    >>> rprod((1 - x, 1 - x**2, 1 - x**5), 8)
    x**7 + x**6 - x**5 + x**3 - x**2 - x + 1
    """
    from sympy.polys.ring_series import rs_mul
    if prec is None:
        p = f[0].ring(1)
        for i in f:
            p *= i
        return p
    x = f[0].ring.x
    p = f[0].ring(1)
    for i in f:
        p = rs_mul(p, i, x, prec)
    return p

>>> from sympy import ZZ
>>> from sympy.polys.ring_series import rs_series_inversion
>>> R, x = ring('x', ZZ)
>>> d = rprod([1-x**i for i in (1,5,10,25,50,100)])
>>> rs_series_inversion(d, x, 101)
293*x**100 + 252*x**99 + 252*x**98 + 252*x**97 + 252*x**96 + 252*x**95 + 218*x**94 + 218*x**93 + 218*x**92 + 218*x**91 + 218*x**90 + 187*x**89 + 187*x**88 + 187*x**87 + 187*x**86 + 187*x**85 + 159*x**84 + 159*x**83 + 159*x**82 + 159*x**81 + 159*x**80 + 134*x**79 + 134*x**78 + 134*x**77 + 134*x**76 + 134*x**75 + 112*x**74 + 112*x**73 + 112*x**72 + 112*x**71 + 112*x**70 + 93*x**69 + 93*x**68 + 93*x**67 + 93*x**66 + 93*x**65 + 77*x**64 + 77*x**63 + 77*x**62 + 77*x**61 + 77*x**60 + 62*x**59 + 62*x**58 + 62*x**57 + 62*x**56 + 62*x**55 + 50*x**54 + 50*x**53 + 50*x**52 + 50*x**51 + 50*x**50 + 39*x**49 + 39*x**48 + 39*x**47 + 39*x**46 + 39*x**45 + 31*x**44 + 31*x**43 + 31*x**42 + 31*x**41 + 31*x**40 + 24*x**39 + 24*x**38 + 24*x**37 + 24*x**36 + 24*x**35 + 18*x**34 + 18*x**33 + 18*x**32 + 18*x**31 + 18*x**30 + 13*x**29 + 13*x**28 + 13*x**27 + 13*x**26 + 13*x**25 + 9*x**24 + 9*x**23 + 9*x**22 + 9*x**21 + 9*x**20 + 6*x**19 + 6*x**18 + 6*x**17 + 6*x**16 + 6*x**15 + 4*x**14 + 4*x**13 + 4*x**12 + 4*x**11 + 4*x**10 + 2*x**9 + 2*x**8 + 2*x**7 + 2*x**6 + 2*x**5 + x**4 + x**3 + x**2 + x + 1
>>> _.coeff(x**100)
293

There is an SO question about this.

[smichr]

See also this post by JohanC regarding generating functions though that approach does not work for the particular case of the OP.