How to define a variable function in advance

General Usage
1

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)

2

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
3

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?

4

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
5

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.

6

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
7

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) ?

8

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 GitHub - JuliaAlgebra/StaticPolynomials.jl: Fast evaluation of multivariate polynomials package.

I donā€™t understand your question about g(x) * f(x).

9

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.