diff --git a/Project.toml b/Project.toml
index 59f35223..268eef8b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.12.2"
+version = "0.13.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/intervals.jl b/src/intervals.jl
index 522e20c7..ab3e3fca 100644
--- a/src/intervals.jl
+++ b/src/intervals.jl
@@ -1,15 +1,151 @@
 using .IntervalArithmetic
 
 # Method used for Taylor1{Interval{T}}^n
-function ^(a::Taylor1{T}, n::Integer) where {T<:Interval}
-    n == 0 && return one(a)
-    n == 1 && return copy(a)
-    n == 2 && return square(a)
-    n < 0 && return a^float(n)
-    return power_by_squaring(a, n)
+for T in (:Taylor1, :TaylorN)
+    @eval function ^(a::$T{Interval{S}}, n::Integer) where {S<:Real}
+        n == 0 && return one(a)
+        n == 1 && return copy(a)
+        n == 2 && return square(a)
+        n < 0 && return a^float(n)
+        return power_by_squaring(a, n)
+    end
 end
 
+function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
+    a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
+    aux = one(a0)^r
+
+    iszero(r) && return Taylor1(aux, a.order)
+    aa = one(aux) * a
+    aa[0] = one(aux) * a0
+    r == 1 && return aa
+    r == 2 && return square(aa)
+    r == 1/2 && return sqrt(aa)
+
+    l0 = findfirst(a)
+    lnull = trunc(Int, r*l0 )
+    if (a.order-lnull < 0) || (lnull > a.order)
+        return Taylor1( zero(aux), a.order)
+    end
+    c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
+    c = Taylor1(zero(aux), c_order)
+    for k = 0:c_order
+        pow!(c, aa, r, k)
+    end
+
+    return c
+end
+function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
+    a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
+    a0r = a0^r
+    aux = one(a0r)
+
+    iszero(r) && return TaylorN(aux, a.order)
+    aa = aux * a
+    aa[0] = aux * a0
+    r == 1 && return aa
+    r == 2 && return square(aa)
+    r == 1/2 && return sqrt(aa)
+    isinteger(r) && return aa^round(Int,r)
+
+    # @assert !iszero(a0)
+    iszero(a0) && throw(DomainError(a,
+        """The 0-th order TaylorN coefficient must be non-zero
+        in order to expand `^` around 0."""))
+
+    c = TaylorN( a0r, a.order)
+    for ord in 1:a.order
+        pow!(c, aa, r, ord)
+    end
+
+    return c
+end
+
+
+for T in (:Taylor1, :TaylorN)
+    @eval function log(a::$T{Interval{S}}) where {S<:Real}
+        iszero(constant_term(a)) && throw(DomainError(a,
+                """The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
+        a0 = constant_term(a) ∩ Interval(S(0.0), S(Inf))
+        order = a.order
+        aux = log(a0)
+        aa = one(aux) * a
+        aa[0] = one(aux) * a0
+        c = $T( aux, order )
+        for k in eachindex(a)
+            log!(c, aa, k)
+        end
+        return c
+    end
 
+    @eval function asin(a::$T{Interval{S}}) where {S<:Real}
+        a0 = constant_term(a) ∩ Interval(-one(S), one(S))
+        a0^2 == one(a0) && throw(DomainError(a,
+            """Series expansion of asin(x) diverges at x = ±1."""))
+
+        order = a.order
+        aux = asin(a0)
+        aa = one(aux) * a
+        aa[0] = one(aux) * a0
+        c = $T( aux, order )
+        r = $T( sqrt(1 - a0^2), order )
+        for k in eachindex(a)
+            asin!(c, aa, r, k)
+        end
+        return c
+    end
+
+    @eval function acos(a::$T{Interval{S}}) where {S<:Real}
+        a0 = constant_term(a) ∩ Interval(-one(S), one(S))
+        a0^2 == one(a0) && throw(DomainError(a,
+        """Series expansion of asin(x) diverges at x = ±1."""))
+
+        order = a.order
+        aux = acos(a0)
+        aa = one(aux) * a
+        aa[0] = one(aux) * a0
+        c = $T( aux, order )
+        r = $T( sqrt(1 - a0^2), order )
+        for k in eachindex(a)
+            acos!(c, aa, r, k)
+        end
+        return c
+    end
+
+    @eval function acosh(a::$T{Interval{S}}) where {S<:Real}
+        a0 = constant_term(a) ∩ Interval(one(S), S(Inf))
+        a0^2 == one(a0) && throw(DomainError(a,
+            """Series expansion of acosh(x) diverges at x = ±1."""))
+
+        order = a.order
+        aux = acosh(a0)
+        aa = one(aux) * a
+        aa[0] = one(aux) * a0
+        c = $T( aux, order )
+        r = $T( sqrt(a0^2 - 1), order )
+        for k in eachindex(a)
+            acosh!(c, aa, r, k)
+        end
+        return c
+    end
+
+    @eval function atanh(a::$T{Interval{S}}) where {S<:Real}
+        order = a.order
+        a0 = constant_term(a) ∩ Interval(-one(S), one(S))
+        aux = atanh(a0)
+        aa = one(aux) * a
+        aa[0] = one(aux) * a0
+        c = $T( aux, order)
+        r = $T(one(aux) - a0^2, order)
+        iszero(constant_term(r)) && throw(DomainError(a,
+            """Series expansion of atanh(x) diverges at x = ±1."""))
+
+        for k in eachindex(a)
+            atanh!(c, aa, r, k)
+        end
+        return c
+    end
+end
 
 function evaluate(a::Taylor1, dx::Interval)
     order = a.order
diff --git a/src/power.jl b/src/power.jl
index d50d2ae5..eeccc443 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -85,11 +85,8 @@ end
 
 ## Real power ##
 function ^(a::Taylor1, r::S) where {S<:Real}
-    # println()
     a0 = constant_term(a)
-    # @show(a, a0)
     aux = one(a0)^r
-    # @show(aux)
 
     iszero(r) && return Taylor1(aux, a.order)
     aa = aux*a
@@ -103,12 +100,11 @@ function ^(a::Taylor1, r::S) where {S<:Real}
         return Taylor1( zero(aux), a.order)
     end
     c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
-    # @show(c_order)
     c = Taylor1(zero(aux), c_order)
     for k = 0:c_order
         pow!(c, aa, r, k)
     end
-    # println()
+
     return c
 end
 
diff --git a/test/intervals.jl b/test/intervals.jl
index a340f170..18d2e93c 100644
--- a/test/intervals.jl
+++ b/test/intervals.jl
@@ -99,4 +99,42 @@ eeuler = Base.MathConstants.e
     displayBigO(true)
     @test string(Taylor1(vc, 5)) ==
         " ( [1.5, 2] + [0, 0]im ) - ( [1, 2] + [-1, 1]im ) t + ( [-1, 1.5] + [-1, 1.5]im ) t² + ( [0, 0] + [-1, 1.5]im ) t³ + 𝒪(t⁶)"
+
+    # Tests related to Iss #311
+    # `sqrt` and `pow` defined on Interval(0,Inf)
+    @test_throws DomainError sqrt(ti)
+    @test sqrt(Interval(0.0, 1.e-15) + ti) == sqrt(Interval(-1.e-15, 1.e-15) + ti)
+    aa = sqrt(sqrt(Interval(0.0, 1.e-15) + ti))
+    @test aa == sqrt(sqrt(Interval(-1.e-15, 1.e-15) + ti))
+    bb = (Interval(0.0, 1.e-15) + ti)^(1/4)
+    @test bb == (Interval(-1.e-15, 1.e-15) + ti)^(1/4)
+    @test all(aa.coeffs[2:end] .⊂ bb.coeffs[2:end])
+    @test_throws DomainError sqrt(x)
+    @test sqrt(Interval(-1,1)+x) == sqrt(Interval(0,1)+x)
+    @test (Interval(-1,1)+x)^(1/4) == (Interval(0,1)+x)^(1/4)
+
+    # `log` defined on Interval(0,Inf)
+    @test_throws DomainError log(ti)
+    @test log(Interval(0.0, 1.e-15) + ti) == log(Interval(-1.e-15, 1.e-15) + ti)
+    @test_throws DomainError log(y)
+    @test log(Interval(0.0, 1.e-15) + y) == log(Interval(-1.e-15, 1.e-15) + y)
+    # `asin` and `acos` defined on Interval(-1,1)
+    @test_throws DomainError asin(Interval(1.0 .. 2.0) + ti)
+    @test asin(Interval(-2.0 .. 0.0) + ti) == asin(Interval(-1,0) + ti)
+    @test_throws DomainError acos(Interval(1.0 .. 2.0) + ti)
+    @test acos(Interval(-2.0 .. 0.0) + ti) == acos(Interval(-1,0) + ti)
+    @test_throws DomainError asin(Interval(1.0 .. 2.0) + x)
+    @test asin(Interval(-2.0 .. 0.0) + x) == asin(Interval(-1,0) + x)
+    @test_throws DomainError acos(Interval(1.0 .. 2.0) + x)
+    @test acos(Interval(-2.0 .. 0.0) + x) == acos(Interval(-1,0) + x)
+    # acosh defined on Interval(1,Inf)
+    @test_throws DomainError acosh(Interval(0.0 .. 1.0) + ti)
+    @test acosh(Interval(0.0 .. 2.0) + ti) == acosh(Interval(1.0 .. 2.0) + ti)
+    @test_throws DomainError acosh(Interval(0.0 .. 1.0) + x)
+    @test acosh(Interval(0.0 .. 2.0) + x) == acosh(Interval(1.0 .. 2.0) + x)
+    # atanh defined on Interval(-1,1)
+    @test_throws DomainError atanh(Interval(1.0 .. 1.0) + ti)
+    @test atanh(Interval(-2.0 .. 0.0) + ti) == atanh(Interval(-1.0 .. 0.0) + ti)
+    @test_throws DomainError atanh(Interval(1.0 .. 1.0) + y)
+    @test atanh(Interval(-2.0 .. 0.0) + y) == atanh(Interval(-1.0 .. 0.0) + y)
 end