Ok I found a solution using a branch of `TaylorSeries.jl`

that had almost everything I needed. Take a look at this code :

```
# ] add TaylorSeries#lb/Symbolics
using Symbolics, TaylorSeries, BenchmarkTools
############### Parameters
N = 200
m = 10
p,α,s = (rand(N), rand(N), rand())
p ./= sum(p)
################ Old version
function f(t,p,α,s)
x = (t-1) / (t*(1-s) - (1+s))
rez = zero(x)
for i in eachindex(p)
rez += p[i] * x^α[i]
end
return rez * sqrt(2) / (1-t)
end
const tayl = Taylor1(m)
coefficients(args...) = f(tayl,args...).coeffs
##############"" New version
@variables t₀ αₛ sₛ;
tₛ = t₀ + Taylor1(m);
tₛ = tₛ - t₀;
function fi(t,α,s)
x = (t-1) / (t -t*s - 1-s)
return x^α / (1-t)
end
# This errors as the power is not defined.
# fi(tₛ,αₛ,sₛ)
# Let's define it, according to the definition in TaylorSeries.jl
import Base.^
^(a::Taylor1{T}, r::Num) where {T<:Number} = exp(r * log(a)) # quick dirty workaround.
# Let's now extract coefficients from the symbolic Taylor series, at t₀ = 0
symbolic_coefs = fi(tₛ,αₛ,sₛ).coeffs;
# Let's make a few more substitutions :
symbolic_coefs = substitute.(symbolic_coefs,(Dict(log(-1/(-1-sₛ)) => -log(1-sₛ)),));
symbolic_coefs = substitute.(symbolic_coefs,(Dict(exp(-αₛ*log(1 - sₛ)) => (1/(1-sₛ))^αₛ),));
#Compile
coefs_expr= build_function(symbolic_coefs,αₛ,sₛ);
taylor_fi = eval(coefs_expr[1]);
coefficients2(p,α,s) = sqrt(2) * sum(pj .* taylor_fi(αj,s) for (pj,αj) in zip(p,α))
coefficients(p,α,s)
coefficients2(p,α,s)
@btime coefficients(p,α,s); # 52.500 μs (822 allocations: 114.52 KiB)
@btime coefficients2(p,α,s); # 101.800 μs (1601 allocations: 106.27 KiB)... Epic fail.
```

However, there is still a lot of redundancy is the symbolic expression for the vector `symbolic_coefs`

. Do you know how I could produce a function that exploits this redundancy, or is it just not possible yet ? I did a little by hand through the `substitute`

calls, but this is not enough, clearly.

For example, the expression of `symbolic_coefs`

contains a lot of `1 - sₛ`

, `1 + sₛ`

and fraction of those two values all over. How could I tell `Symbolics.jl`

that, when compiling the function, it should probably pre-compute these two values, their inverses and their ratio to drastically reduce the number of divisions ? Is it even possible, or am I asking too much ?