Power ^ allocates

[saschatimme]

It seems that x^2 currently allocates, which I think should not happen since Base has a special handling for these cases.

julia> using IntervalArithmetic, BenchmarkTools
julia> x = interval(rand())
[0.375296, 0.375297]
julia> pow(a) = a^2
pow (generic function with 1 method)
julia> @benchmark pow($x)
BenchmarkTools.Trial: 
  memory estimate:  2.15 KiB
  allocs estimate:  55
  --------------
  minimum time:     4.419 μs (0.00% GC)
  median time:      4.777 μs (0.00% GC)
  mean time:        6.551 μs (20.78% GC)
  maximum time:     2.424 ms (99.35% GC)
  --------------
  samples:          10000
  evals/sample:     7

I found this comment

# overwrite new behaviour for small integer powers:
# ^{p}(x::IntervalArithmetic.Interval, ::Type{Val{p}}) = x^p

So is this just not implemented or is there a specific reason to not include these?

[dpsanders]

This is an unfortunate corner of the package at the moment:

The allocation is due to the fact that the only way to get correct rounding for powers is using MPFR, so this is in fact converting the interval to a BigFloat interval, doing the power, and then converting back :/

However, if you are happy to have a slightly wider interval, there is actually an exported function pow that uses a fast path (power by squaring for integer powers):

julia> @btime ($x)^2
  3.266 μs (51 allocations: 2.02 KiB)
Interval(0.3183098861837905, 0.5000000000000003)

julia> @btime pow($x, 2)
  69.124 ns (0 allocations: 0 bytes)
Interval(0.3183098861837905, 0.5000000000000003)

We had a discussion in #14 about whether to make pow the default behaviour.

In any case, there should certainly be an easy way to choose to use fast powers:
#56

[saschatimme]

The problem with a special function like pow is that I cannot use this library as a drop-in replacement. I would like to provide a fast (default) implementation but also a verified implementation using interval arithmetic without any changes to the code base.
At least I write x^2 only instead of x * x since I know that the former is lowered to the latter. So maybe one could argue that here you should also use the small integer power behaviour from base.

But I understand if you want to resolve this in a more general way. I made a suggestion in #56 for this.

[dpsanders]

Actually you can just redefine ^ in your code before using it:

julia> import Base: ^

julia> ^(x::Interval{Float64}, p::Integer) = pow(x, p)

and similarly for Float64 powers etc. if you need them.

By the way, x^2 is not equal to x*x for intervals.

I think the small power business has gone away in master?

[rdeits]

This came up on Discourse with a pretty significant performance impact: https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/3?u=rdeits

[rdeits]

By the way, the special handling of literal powers is still present in Julia 1.0, so it's totally possible to intercept integer literal powers and dispatch them to pow. This could, I think, just be done inside IntervalArithmetic.jl:

julia> using IntervalArithmetic, BenchmarkTools

julia> x = Interval(1.0, 1.0)
[1, 1]

julia> @btime $x^2
  1.327 μs (51 allocations: 2.03 KiB)
[1, 1]

# This definition could live in IntervalArithmetic.jl
julia> Base.literal_pow(::typeof(^), x::Interval, ::Val{N}) where {N} = pow(x, N)

julia> @btime $x^2
  37.049 ns (0 allocations: 0 bytes)
[1, 1]

[stevengj]

You could also intercept small literal powers and dispatch them to an optimal sequence of inlined multiplications.

[dpsanders]

This is definitely a bad situation. I think we should now reverse pow and ^ so that ^ is fast and pow is optimally accurate.

Cc @lbenet

[lbenet]

At least, it is worth giving it a try.

[lucaferranti]

should be fixed in #388

[OlivierHnt]

Solved by PR #593.