# How to define a variable function in advance

I want to define families of polynomial functions for which I can compute the coefficients in advance. For example, (here I stop at degree 3)

``````function poly1(n::Int64, x)
if n == 0
return 1.
elseif n == 1
return x
elseif n == 2
return x^2 +0.595238x-7.875
elseif n == 3
return x^3 + 1.62521x^2 -12.8858x -8.111
end
end
``````

I can also do a more generic code like

``````function poly2(n::Int64, x,  coeff)
if n == 0
return 1.
else
X = [x^j for j in 0:n]
return coeff[n]'*X
end
end
``````

where `coeff` is an array of array here
`coeff = [[0, 1.0],[-7.875, 0.595238, 1.0],[-8.111, -12.8858, 1.62521, 1.0] ]`

The first one is far more efficient

``````using BenchmarkTools
@btime poly1(0,1.3);
0.021 ns (0 allocations: 0 bytes)
@btime poly1(1,1.3);
0.021 ns (0 allocations: 0 bytes)
@btime poly1(2,1.3);
0.021 ns (0 allocations: 0 bytes)

@btime poly2(0,1.3,coeff);
26.810 ns (1 allocation: 16 bytes)
@btime poly2(1,1.3,coeff);
94.884 ns (2 allocations: 112 bytes)
@btime poly2(2,1.3,coeff);
103.483 ns (2 allocations: 128 bytes)
``````

My question is how to define automatically in advance (once I know the `coeff` array), my polynomial functions for `n=0:nmax` as in `poly1` the most efficiently? (Because I need to call these functions a lot of time after)

Note that this is a terrible way to evaluate polynomials. Use `@evalpoly(x, -8.111, -12.8858, 1.62521, 1)`, which basically expands to Hornerās rule.

``````X = [x^j for j in 0:n]
return coeff[n]'*X
``````

This is also a terrible way to evaluate polynomials, even given a coefficient array. Use Hornerās rule. (If you use the Polynomials.jl package to store polynomials, it will do this. Or just use `evalpoly` in Julia 1.4.)

You can use code generation to generate functions for given polynomials, but I would do this as a last resort. First, fix your code.

In fact, if you store the polynomial coefficients as tuples rather than arrays, the `evalpoly` function, which I believe was added in Julia 1.4, should be basically as fast as inlined code generation.

6 Likes

I donāt have Julia v1.4 (still on v1.0). I can do

``````function poly1_correct(n::Int64, x)
if n == 0
return 1.
elseif n == 1
return x
elseif n == 2
return @evalpoly(x,-7.875,0.595238,1)
elseif n == 3
return @evalpoly(x, -8.111,-12.8858,1.62521,1)
end
end
``````

For the one with `coeff` I am not sur how to code that
I tried

``````for n in 1:3
@eval begin
function g(n,x)
@evalpoly x \$(coeff[n]...)
end
end
end
``````

it works fast
`0.021 ns (0 allocations: 0 bytes)`
but it does not work because it returns only `g(3,x)`. Note that I could not get it working without the first `@eval`.
Maybe with Julia 1.4 it will work more smoothly?

You can always implement Hornerās method yourselfā¦

Iām not sure why you have to use Julia 1.0 and not 1.4, but you can easily adapt the implementation of `evalpoly` from Julia 1.4 into your Julia 1.0 code: https://github.com/JuliaLang/julia/blob/712731e1a1e9e85a272343ba2e773707ca90a466/base/math.jl#L101-L193

For example, the simplest generic method is:

``````# evaluate p(x), where the p argument is the coefficients of a polynomial
function evalpoly(x, p)
N = length(p)
ex = p[end]
for i in N-1:-1:1
ex = muladd(x, ex, p[i])
end
ex
end
``````
3 Likes

Thanks for the tip. I agree that Hornerās rule is better:

``````function poly_3(n::Int64, x, coeff)
if n == 0
return 1.
else
return evalpoly(x,coeff[n])
end
end

@btime poly_3(0,1.3,coeff);
@btime poly_3(1,1.3,coeff);
@btime poly_3(2,1.3,coeff);
@btime poly_3(3,1.3,coeff);
27.258 ns (1 allocation: 16 bytes)
34.941 ns (1 allocation: 16 bytes)
37.842 ns (1 allocation: 16 bytes)
37.923 ns (1 allocation: 16 bytes)
``````

But it still does not match the timing of `poly1` function defined at the top of the thread.
Is there a way to define these functions coefficients before usage?
It is not only for this problem but a more generic questions (I have polynomial product with fixed coefficients that I also want to pre-evaluate because I need to call them a lot of time).
You mentioned code generation I tried to look it up but I am not sure how to apply it here, any advices?
Thanks for the help.

Yes, because that inlines/unrolls the loop. Thatās the advantage of using tuples (whose size is known to the compiler) with the new `evalpoly` implementation in Julia 1.4.

2 Likes

I installed Julia 1.4

``````@btime evalpoly(1.3,(-8.111, -12.8858, 1.62521, 1.0))
0.021 ns (0 allocations: 0 bytes)
``````

but

`````` c3 = (-8.111, -12.8858, 1.62521, 1.0)
@btime evalpoly(1.3,c3)
22.939 ns (1 allocation: 16 bytes)
``````

How are we suppose to define these tuples without writting all the coefficients explicitly? I cannot find an example of that.

Also, my other question still holds, how when you have two (or more) functions `f(x)` and `g(x)` can you define and āpre-compileā `h(x) = g(x)*f(x)` ?

You are not benchmarking correctly; you need to āinterpolateā to avoid accessing the global variable `c3`. Unfortunately currently you need to do this as follows with the strange `Ref` syntax and the `[]` indexing:

``````julia> @btime evalpoly(1.3, \$(Ref(c3))[])
1.652 ns (0 allocations: 0 bytes)
-19.9189351
``````

The `0.021 ns` in the first benchmark suggests / indicates that the compiler is optimizing the whole calculation to just calculating the result.

Where do you want to take the coefficients from?
You may want to check the https://github.com/JuliaAlgebra/StaticPolynomials.jl package.

I donāt understand your question about `g(x) * f(x)`.

If you search the forums on the document about performance and global variables you will find descriptions of why accessing global variables has bad performance and how to test around it. A simple way to do it is:

``````@btime evalpoly(1.3, \$c3)
``````

However like @dpsanders said, if the time is less than a nano second the compiler is doing optimizations on you and you are not benchmarking the actual operation.