It’s still not clear to me what you are looking for. You say

GetCoeff, that takes an integer n, a power series z and a monomial x^n,

But, the function you posted `GetCoeff`

takes two arguments, not three; and they both appear to be integers, rather than a power series or monomial. You refer once to `x`

; in `x^n`

. You refer to a fixed polynomial z + z^2 + z^3. The function `GetCoeff`

returns a polynomial in `z[i]`

for i \in [0,1,...n]. Do these refer to the same z ? Perhaps `i`

represents integer-valued t ? It would help if you explain things more precisely.

In general, these algebra tools involve a trade-off between flexibility-expressiveness on one hand and efficiency on the other. As you noted, Mathematica is more expressive than efficient. (Although a few decades and a few hundred million dollars can go a long way to closing this gap.) I’m not sure, but it looks like Nemo.jl and AbstractAlgebra.jl can’t do what you want easily; they want numeric coefficients. But, maybe asking the developers is worthwhile. There is also JuliaAlgebra, as well as a few other polynomial-related packages.

@dpsanders suggested TaylorSeries.jl, which looks like a more likely candidate than the others.

Following is your Mathematica code translated into Symata.jl. But, I don’t have access to Mathematica at the moment, so I can’t test it.

*getcoeff.sj*

```
GetCoeff(p_, n_) := GetCoeff(p, n) = Expand(If(
p == 1, z(n),
If(p == 2, Sum(z(i) * z(n - i), [i, 0, n]),
Apply(Plus, Table(GetCoeff(p - 1, i) * z(n - i), [i, 0, n])))))
```

```
symata 12> Get("getcoeff.sj")
symata 12> GetCoeff(3,5)
Out(12) = 7z(0)*z(2)*z(3) + 7z(0)*z(1)*z(4) + 7*z(0)^2*z(5)
symata 13> ? GetCoeff
GetCoeff(2,0)=z(0)^2
GetCoeff(2,1)=(2z(0)*z(1))
GetCoeff(2,2)=(3z(0)*z(2))
GetCoeff(2,3)=(4z(0)*z(3))
GetCoeff(2,4)=(5z(0)*z(4))
GetCoeff(2,5)=(6z(0)*z(5))
GetCoeff(3,5)=(7z(0)*z(2)*z(3) + 7z(0)*z(1)*z(4) + 7*z(0)^2*z(5))
GetCoeff(p_,n_):=GetCoeff(p,n)=CompoundExpression(nothing,Expand(If(p == 1,z(n),If(p == 2,Sum((z(i)*z(n - i)),[i,0,n]),Apply(Plus,Table((GetCoeff((p - 1),i)*z(n - i)),[i,0,n]))))))
```

You can probably gain some efficiency by using `Symata`

from the `Julia`

side so that you can control the evaluation. Or, you can work directly in SymPy.jl. Or maybe SymEngine.jl, but it is not well-developed compared to `SymPy`

.