Running a search algorithm

edit

You can speed up things by factoring the denominator of $f(x,y)$ and see if for each factor you can find its multiple of the form $1+x^l+y^m$ - these will give you candidate denominators in the decomposition. Given that $1+x^l+y^m$ is irreducible (over $\mathbb{Q}$ at least), each factor must simply be of this form.

A "brute-force" approach to this problem is doubtful :

Such an approach would generate all possible sum of 16 terms, each representing one quadruplet $i, j, k, l \in \{0, \dots, 5\}$. Since there are $6^4=1296$ such potential summation terms, there are ${6^4}^{16}$ possible 16-terms sums.But :

sage: ((6^4)^16).log(10).n()
49.8016800245532

meaning that the number of sums is about $6.44\cdot 10^{49}$, which seems a bit too large for systematic exploration...

Limiting the size of the exploration by other means is necessary. But this is the feathers of a totally different horse...

HTH, but doubting it,

Prehaps coud itertools.product be of some help:

from itertools import product
P = product(range(6), repeat=4)
for i,j,k,l in P:
    print(i,j,k,l)