Expanding a bivariate exponential generating function

edit

Expanding an univariate exponential generating function can be done like this:

def egfExpand1(f, size):
    x = var('x')
    return taylor(f(x), x, 0, size).power_series(SR).egf_to_ogf().list()

For example egfExpand1(sec, 10) returns [1, 0, 1, 0, 5, 0, 61, 0, 1385, 0, 50521].

But how can I expand a bivariate exponential generating function? Say

def f(x, y): return exp(x * y) * sec(x)

def egfExpand2(f, size):
    return ...

The expected output is an integer triangle (i.e. a list of integer lists).

The example would return an unsigned version of A119879, which starts:

                           1
                          0, 1
                        1, 0, 1
                       0, 3, 0, 1
                     5, 0, 6, 0, 1

Edit:

Frédéric suggested the following solution, slightly rewritten here.

def egfExpand2(f, size):
    y = polygen(QQ, "y")
    x = LazyPowerSeriesRing(y.parent(), "x").gen()
    return [list(f(x, y)[n] * factorial(n)) for n in range(size)]

f = lambda x, y: exp(x * y) * sec(x)
egfExpand2(f, 10)

The univariate case can also be written more elegantly with this method:

def egfExpand1(f, size: int):
    x = LazyPowerSeriesRing(QQ, "x").gen()
    return [f(x)[n] * factorial(n) for n in range(size)]

egfExpand1(sec, 11)