Ok, I’m not sure whether anyone is actually interested in this kind of problems, but I thought I should post my solution just in case. The code is rather hacky but works better than what I expected.

I ended up using SymEngine.jl because I don’t know how else to produce the equations.

It calculates the first 100 coefficients in about 0.14s and the first 40 in less than 10s! That may not sound too impressive, but you need to keep in mind that the coefficients grow as factorial so the numerators and denominators are ridiculously large numbers. Anyway, here’s the code:

```
using SymEngine
using Memoize
mutable struct InductivePowerSeries
coeffs::Array{Rational{BigInt},1}
startingOrder::Int64
knownUpTo::Int64
symbol::SymEngine.Basic
dummy::SymEngine.Basic
end
function InductivePowerSeries(N::Int,s::String)
eval(parse(s*" = symbols(\""*s*"\")"))
eval(parse(s*"0 = symbols(\""*s*"0\")"))
return InductivePowerSeries(Array{Rational{BigInt},1}(N),1,0,eval(s),eval(s*"0"))
end
function getCoeff(s::InductivePowerSeries,n::Int)::Union{Rational{BigInt},SymEngine.Basic}
if n < s.startingOrder
return Rational{BigInt}(0)
end
if n <= s.knownUpTo
return s.coeffs[n]
end
if n == s.knownUpTo + 1
return s.symbol
end
if n == s.knownUpTo + 2
return s.dummy
end
error("Error: the power series is known up to order $(s.knownUpTo) so the term at position $n cannot be returned.")
end
function setNextCoeff!(s::InductivePowerSeries,n::Int,value::Rational{BigInt})
if n == s.knownUpTo + 1
s.coeffs[n] = value
s.knownUpTo += 1
return
else
error("Error: the power series is known up to order $(s.knownUpTo) so the term at position $n cannot be set.")
end
end
setNextCoeff!(s::InductivePowerSeries,n::Int,value::Union{Int,Rational{Int},BigInt}) = setNextCoeff!(s,n,Rational{BigInt}(value))
function getCoeffOfPower(s::InductivePowerSeries,p::Int,n::Int)::Union{Rational{BigInt},SymEngine.Basic}
if n < p * s.startingOrder
return Rational{BigInt}(0)
end
if n < p * s.startingOrder + s.knownUpTo
return getCoeffOfPowerNumber(s,p,n)
end
if p == 1
return getCoeff(s,n)
end
if n <= p + s.knownUpTo + 1
sum = 0*s.symbol
for i = s.startingOrder*(p-1):n-s.startingOrder
sum += getCoeffOfPower(s,p-1,i)*getCoeff(s,n - i)
end
return sum
end
error("Error: the power series is known up to order $(s.knownUpTo) and the order $n of power $p was asked.")
end
@memoize function getCoeffOfPowerNumber(s::InductivePowerSeries,p::Int,n::Int)::Rational{BigInt}
if p == 1
return getCoeff(s,n)
elseif p == 2
sum = Rational{BigInt}(0)
for i = s.startingOrder:n-s.startingOrder
sum += getCoeff(s,i)*getCoeff(s,n-i)
end
return sum
else
sum = Rational{BigInt}(0)
for i = s.startingOrder*(p-1):n-s.startingOrder
sum += getCoeffOfPowerNumber(s,p-1,i)*getCoeff(s,n-i)
end
return sum
end
end
function translationCoeff(p::Int,n::Int)::BigInt
if (p + n) % 2 == 0
return binomial(BigInt(n-1),p-1)
else
return -binomial(BigInt(n-1),p-1)
end
end
function getTranslatedCoeffNumber(s::InductivePowerSeries,n::Int)::Rational{BigInt}
sum = BigInt(0)
for i = 1:n
sum += translationCoeff(i,n)*getCoeff(s,i)
end
return sum
end
function getTranslatedCoeffButTheLastNumber(s::InductivePowerSeries,n::Int)::Rational{BigInt}
sum = BigInt(0)
for i = 1:n-1
sum += translationCoeff(i,n)*getCoeff(s,i)
end
return sum
end
function getTranslatedCoeffButTheLast2Number(s::InductivePowerSeries,n::Int)::Rational{BigInt}
sum = BigInt(0)
for i = 1:n-2
sum += translationCoeff(i,n)*getCoeff(s,i)
end
return sum
end
function getTranslatedCoeff(s::InductivePowerSeries,n::Int)::Union{Rational{BigInt},SymEngine.Basic}
if n > s.knownUpTo + 2
error("Error: the power series is known up to order $(s.knownUpTo) but the translation coefficient asked for is at order $n")
end
if n <= s.knownUpTo
return getTranslatedCoeffNumber(s,n)
end
if n == s.knownUpTo + 1
return s.symbol + getTranslatedCoeffButTheLastNumber(s,n)
end
if n == s.knownUpTo + 2
return s.dummy-(n-1)*s.symbol + getTranslatedCoeffButTheLast2Number(s,n)
end
end
function calculateCoeffs(Z::InductivePowerSeries,MaxDeg::Int)
for n = Z.knownUpTo+2:MaxDeg + 1
p = expand(getTranslatedCoeff(Z,n) - getCoeffOfPower(Z,1,n) - getCoeffOfPower(Z,2,n) - getCoeffOfPower(Z,3,n))
value = -Rational{BigInt}(p(0))//Rational{BigInt}((p-p(0))(1))
setNextCoeff!(Z,n-1,value)
end
end
```

In case someone really cares to understand the code, I can try to explain it, but I would bet otherwise. Thank you all for the discussion, it was actually very useful.