# 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
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.