Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add pow for interval^real and interval^interval #42

Merged
merged 3 commits into from
Jun 11, 2017
Merged

Add pow for interval^real and interval^interval #42

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c6e38fa
Add pow for interval^real and interval^interval
dpsanders May 27, 2017
8aed052
Add some tests for real and integer powers using pow
dpsanders Jun 7, 2017
40eb21a
Fix subset
dpsanders Jun 7, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 13 additions & 5 deletions src/intervals/functions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,18 +191,18 @@ function sqrt{T}(a::Interval{T})
@round(sqrt(a.lo), sqrt(a.hi)) # `sqrt` is correctly-rounded
end

doc"""
"""
pow(x::Interval, n::Integer)

A faster implementation of `x^n` using `power_by_squaring`.
`pow(x, n) will usually return an interval that is slightly larger than that calculated by `x^n`, but is guaranteed to be a correct
A faster implementation of ``x^n``, currently using `power_by_squaring`.
`pow(x, n)` will usually return an interval that is slightly larger than that calculated by `x^n`, but is guaranteed to be a correct
enclosure when using multiplication with correct rounding.
"""
function pow{T}(x::Interval{T}, n::Integer) # fast integer power
function pow(x::Interval, n::Integer) # fast integer power

isempty(x) && return x

if iseven(n) && zero(T) ∈ x
if iseven(n) && 0 ∈ x

return hull(zero(x),
hull(Base.power_by_squaring(Interval(mig(x)), n),
Expand All @@ -218,6 +218,14 @@ function pow{T}(x::Interval{T}, n::Integer) # fast integer power

end

function pow(x::Interval, y::Real) # fast real power, including for y an Interval

isempty(x) && return x

return exp(y * log(x))

end




Expand Down
72 changes: 51 additions & 21 deletions test/interval_tests/numeric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,33 +194,63 @@ end
end

@testset "Fast power" begin
x = 1..2
@test pow(x, 2) == pow(-x, 2) == Interval(1, 4)
@test pow(-x, 3) == Interval(-8.0, -1.0)

@test pow(-1..2, 2) == 0..4
@test pow(-1..2, 3) == -1..8
@test pow(-1..2, 4) == 0..16
@testset "Fast integer powers" begin
x = 1..2
@test pow(x, 2) == pow(-x, 2) == Interval(1, 4)
@test pow(-x, 3) == Interval(-8.0, -1.0)

@test pow(@biginterval(-1, 2), 2) == 0..4
@test pow(@biginterval(-1, 2), 3) == -1..8
@test pow(@biginterval(1, 2), 2) == 1..4
@test pow(-1..2, 2) == 0..4
@test pow(-1..2, 3) == -1..8
@test pow(-1..2, 4) == 0..16

@test pow(@biginterval(-1, 2), 2) == 0..4
@test pow(@biginterval(-1, 2), 3) == -1..8
@test pow(@biginterval(1, 2), 2) == 1..4

x = @interval(pi)
@test x^100 ⊆ pow(x, 100)
@test x^50 ⊆ pow(x, 50)
@test isinterior(x^50, pow(x, 50))
end
x = @interval(pi)
@test x^100 ⊆ pow(x, 100)
@test x^50 ⊆ pow(x, 50)
@test isinterior(x^50, pow(x, 50))

x = Interval(2)
@test pow(x, 2000) == Interval(realmax(), Inf)
end

@testset "Fast real powers" begin
x = 1..2
@test pow(x, 0.5) == Interval(1.0, 1.4142135623730951)
@test pow(x, 0.5) == x^0.5

@test pow(x, 17.3) != x^17.3

y = 2..3
@test pow(y, -0.5) == Interval(0.5773502691896257, 0.7071067811865476)

y = -2..3
@test pow(y, 2.1) == Interval(0.0, 10.045108566305146)
@test pow(y, 2.1) != y^2.1
@test y^2.1 ⊆ pow(y, 2.1)
end

@testset "Fast interval powers" begin
x = 1..2
@test x^Interval(-1.5, 2.5) == Interval(0.35355339059327373, 5.656854249492381)

@testset "Behaviour near infinity" begin
a = Interval(1e300)
b = Interval(1e9)
y = -2..3
@test pow(y, 2.1) == Interval(0.0, 10.045108566305146)
@test pow(y, Interval(-2, 3)) == Interval(0, Inf)

@test a * b == Interval(realmax(), Inf)
@test pow(y, @interval(2.1)) == Interval(0.0, 10.045108566305146)
end

# comment out for now:
# a = Interval{Float32}(1e38)
# b = Interval{Float32}(1e2)
# @test a * b == Interval{Float32}(realmax(Float32), Inf)
#=
@testset "Float32 intervals" begin

a = Interval{Float32}(1e38)
b = Interval{Float32}(1e2)
@test a * b == Interval{Float32}(realmax(Float32), Inf)
end
=#
end