Powering a Math constant in Julia is slow

Question

I'm pretty sure this has to be a bug...

Ok, so basically, when I pass a Math Constant type to the power (^) function and iterate it several times... The loop is really slow and uses a lot of memory allocation. The example is trivial as it doesnt produce any result, but I had real application where I had this issue.

function test1(N)
    for i = 1:N
        x = pi^4
    end
end

function test2(N)
    for i = 1:N
        x = pi*pi*pi*pi
    end
end

function test3(N)
    pif = float64(pi)
    for i = 1:N
        x = pif^4
    end
end

And this is the result of calling the functions with a big N:

@time test1(10000000)
@time test2(10000000)
@time test3(10000000)

elapsed time: 0.958278949 seconds (320000160 bytes allocated, 56.44% gc time)
elapsed time: 6.341e-6 seconds (80 bytes allocated)
elapsed time: 4.982e-6 seconds (80 bytes allocated)

Is this is a bug? Or is there a logical explanation for it? What would the proper way to power pi to the 4th inside a big loop?

Thanks.

Answer 1

Looking at Constants.jl in Julia 0.3 and 0.4

Although it looks like you are using Julia 0.3, and these results are principally about Julia 0.4, it also has the same problem and the cause is the same which is unnecessary allocation of intermediates.

To get an idea of what is going on, you can type methods(^) and it will produce a list of all of the methods defined for implementation of ^. There are a lot of them so I won't list them, but the important thing is that the Julia file that implements them is listed, and in this case, it's Constants.jl.

Julia 0.3 vs. Julia 0.4

First time in Julia 0.3

julia> @time test1(10000000)
elapsed time: 0.313772825 seconds (320393016 bytes allocated, 37.16% gc time)

In Julia 0.4, owing to the better garbage collector, the problem is diminished, but still present.

julia> @time test1(10000000)
 170.445 milliseconds (20000 k allocations: 305 MB, 6.94% gc time)

julia> @time test2(10000000)
   2.355 microseconds (4 allocations: 144 bytes)

When there is an abnormal number of allocations, always suspect type instability or temporaries.

definition of ^ for pi

In particular, lines 70 through 72 define a generic method for raising a MathConst to a power.

for op in Symbol[:+, :-, :*, :/, :^]
    @eval $op(x::MathConst, y::MathConst) = $op(Float64(x),Float64(y))
end

using ^ for e instead of pi

Note also that later on is a specialization for the constant e, and a new test where pi is replaced by e does not have this same problem.

julia> function test4(N)
           for i = 1:N
               x = e^4
           end
       end
test4 (generic function with 1 method)

run test4

julia> @time test4(10000000)
 108.401 milliseconds (4 allocations: 144 bytes)

conversion of Int to Float64

If a Float is being created from an Int, then this will avoid the problem by starting with a Float.

julia> function test5(N)
           for i = 1:N
               x = pi^4.0
           end
       end

which is what happens

julia> @time test5(10000000)
  65.430 milliseconds (4 allocations: 144 bytes)

define a new set of functions

Finally, one can create a new set of functions (testing in this case), that specifically define Int as the 2nd parameter to avoid the cast. This is not the best approach because of the ambiguity it introduces, but it's good for a test.

julia> for op in Symbol[:+, :-, :*, :/, :^]
           @eval $op(x::MathConst, y::Int64) = $op(Float64(x),y)
       end

Warning: New definition 
    ^(MathConst{sym}, Int64) at none:2
is ambiguous with: 
    ^(MathConst{:e}, Integer) at constants.jl:122.
To fix, define 
    ^(MathConst{:e}, Int64)
before the new definition.

and then re-create test1

julia> function test1(N)
           for i = 1:N
               x = pi^4
           end
       end
test1 (generic function with 1 method)

and then re-run

julia> @time test1(10000000)
   2.757 microseconds (4 allocations: 144 bytes)

do the same thing for Julia 0.3

Same fix, but use float64 instead of Float64

julia> for op in Symbol[:+, :-, :*, :/, :^]
          @eval $op(x::MathConst, y::Int64) = $op(float64(x),y)
       end

and finally re-define and re-run test in in Julia 0.3

julia> @time test1(10000000)
elapsed time: 3.023e-6 seconds (80 bytes allocated)