diff --git a/src/intervals/functions.jl b/src/intervals/functions.jl
index 8fa052518..14cd8960b 100644
--- a/src/intervals/functions.jl
+++ b/src/intervals/functions.jl
@@ -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),
@@ -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
+
 
 
 
diff --git a/test/interval_tests/numeric.jl b/test/interval_tests/numeric.jl
index d7cc39155..178952bd7 100644
--- a/test/interval_tests/numeric.jl
+++ b/test/interval_tests/numeric.jl
@@ -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