Extract coefficients of series of polynomials I have a function $f(x,y)$ which can be formally written as a series of polynomials
$$
f(x,y) = \sum_{n,m}a_{n,m,k}\,p_1(x,y)^n\,p_2(x,y)^m\,p_3(x,y)^k
$$
where 
$$
p_1(x,y) = x^2 + y^2\\
p_2(x,y) = x + y\\
p_3(x,y) = xy\\
$$
The coefficients $a_{n,m,k}$ are unknown and I want to determine them by operating on $f(x,y)$. How can I do that?
For example, the following polynomial
$$
f(x,y)=3 + 3 x + 6 x^2 + 3 x^3 + 3 x^4 + 3 y + 6 x y + 3 x^2 y + 6 x^3 y + 
 6 y^2 + 3 x y^2 + 6 x^2 y^2 + 3 y^3 + 6 x y^3 + 3 y^4
$$
can be written as a sum of polynomials $p_1,p_2,p_3$. How can I derive the coefficients $a_{n,m,k}$?
Edited after Ronald's comment
I am interested in deriving an algorithm to extract the coefficients $a_{n,m,k}$ without expanding the sum up the maximal order of the polynomial and then matching the coefficients. For example, for a function $f(x)$ which admits an expansion around $x=0$, we have
$$
f(x) = \sum_n a_n x^n
$$
and the coefficient $a_n$ can be derived just by acting with derivatives
$$
a_n = \frac{1}{n!} \frac{d^n f(x)}{dx^n}
$$
Analogously, for any function $f(x) = a + b x + c x^2$ the coefficient $c$ can be also derived by looking at the residue at $x=+\infty$
$$
c = -\text{Res}_{x=+\infty}\left[\frac{f(x)}{x^3}\right]
$$
I would like to do something similar. I am interesting if I can derive a formula to extract the coefficients $a_{n,m,k}$ in an analogously way.
 A: The first thing is to check whether the polynomials $p_i(x,y)$ are independent. That is here not the case, i.e., $p_3 = (p_2^2 - p_1)/2$, hence the expansion can be reduced to the simpler form
$$
\sum_{m,n} a_{m,n} (x+y)^m (x^2+y^2)^n
$$
Since the maximum degree of the function $f(x,y)$ is 4, the summation limits can be reduced to
$$
\sum_{n=0}^2 \sum_{m=0}^{4-2 n} a_{m,n} (x+y)^m (x^2+y^2)^n
$$
with at most 9 coefficients $a_{m,n}$. You could now of course expand this in $x$ and $y$ and subtract the function $f(x,y)$. Since this should be zero, you get an linear set of equations for each power $x^i y^j$ in the unknowns $a_{m,n}$, which you can solve.
Note that such a solution might not always exist, i.e., for $f(x,y)=x$ there is no solution.
An alternative but identical approach, is to take all terms of the same order in $f(x,y)$ and make all possible combinations of $p_1$ and $p_2$ that result in the same order. For example for order 3, you find that the function $f(x,y)$ has the following terms
$$
f^{(3)}(x,y) = 3 x^3 + 3 x^2 y + 3 x y^2 + 3 y^3
$$
In the expansion there are only two contributions, i.e., $p_1 p_2$ and $p_2^3$. This can easily be solved by realising that $f^{(3)}(x,y)=(x+y)(x^2+y^2)$.
In a similar fashion you can treat the other orders present in $f(x,y)$ and find the correct expansion.
In view of your clarification, you will have in the more general expression in 2 variables $x$ and $y$
$$
f(x,y) = \sum_{m \geq 0} \sum_{n \geq 0} c_{m,n} x^m y^n
$$
that needs to be matched with an expansion of the type
$$
f(x,y) = \sum_{m \geq 0} \sum_{n \geq 0} a_{m,n} p_1(x,y)^m p_2(x,y)^n
$$
Assuming this can be done (see remark before), you can obtain the coefficient $c_{m,n}$ as 
$$
c_{m,n} = \frac{1}{m! n!} \frac{\partial^m}{\partial x^m}  \frac{\partial^n}{\partial y^n} f(x,y)
$$
By applying this to the second expansion you immediately get a linear equation in the coefficients $a_{m,n}$. In general they will contain several different $a's$ and hence you need to find additional equations to solve them, i.e., they come in groups with $m+n = $ constant (the polynomial order of the term) that need to be solved simultaneously. Note that this is identical with the procedure sketched above.
The generalisation to more variables is straight forward, as well as including more polynomials $p_i$. The later, however, need to be independent in order to give a unique solution.
In principle you could follow another approach, by realising that for $x>y$ you have
$$
x = \frac{p_2}{2} + \frac{\sqrt{2 p_2^2 - p_1}}{2}
$$
$$
y = \frac{p_2}{2} - \frac{\sqrt{2 p_2^2 - p_1}}{2}
$$
Substituting these expressions in $f(x,y)$ will give the required expansion directly.
