multivariate power series computation with recognition

edit

First, a personal comment:

For technical implementation of power series, there is a difference between "lazy" power series and "truncated" power series. It sounds like you're interested in the former, but Sage implements the latter. There are different tradeoffs for the two approaches; please don't suggest that one is "improper" unless you're prepared to make a very thorough case that it is uniformly inferior. To have real weight, such a case would probably need to be accompanied by an actual implementation of the other approach.

Now on to your specific question:

I think the calculations you're interested in can be done just fine with truncated power series. The essential paradigm change is that you need to decide at the outset how many coefficients you're interested in actually computing. This will be the "precision" of the power series that you work with. For multivariable power series, the precision limits the total degree of terms. So in a power series with precision 5, terms like x^3*y^3 would be dropped because they have total degree greater than 5. Some elements can have infinite precision, but only those which are actually polynomials (that is, they have only finitely many nonzero coefficients).

Every power series ring (in Sage) has a default precision -- the precision that it assigns by default to elements that don't already have one. This is only important in your case because exp uses the default precision for the parent ring of its arguments.

Here is some code that would carry out your calculation through total degree 4:

sage: prec = 5                                                                 
sage: R.<x,y> = PowerSeriesRing(QQ, default_prec=prec); R
Multivariate Power Series Ring in x, y over Rational Field
sage: h = sum(factorial(m)*4^n*x^m*y^n for m in range(prec) for n in range(pre>
sage: k = sum((m+n)*x^m*y^n for m in range(prec) for n in range(prec))
sage: exp(x*y)
1 + x*y + 1/2*x^2*y^2 + O(x, y)^5
sage: H = h(x,exp(x*y)*y)
sage: H
1 + x + 4*y + 2*x^2 + 4*x*y + 16*y^2 + 6*x^3 + 8*x^2*y + 20*x*y^2 + 64*y^3 + 24*x^4 + 24*x^3*y + 36*x^2*y^2 + 96*x*y^3 + 256*y^4 + O(x, y)^5
sage: H*k
x + y + 3*x^2 + 7*x*y + 6*y^2 + 7*x^3 + 19*x^2*y + 33*x*y^2 + 27*y^3 + 17*x^4 + 45*x^3*y + 91*x^2*y^2 + 143*x*y^3 + 112*y^4 + 46*x^5 + 119*x^4*y + 219*x^3*y^2 + 407*x^2*y^3 + 605*x*y^4 + 448*y^5 + O(x, y)^6

If you really only need exp(x*y), you could define this just as h and k are defined, and thereby avoid ...