diff --git a/Project.toml b/Project.toml
index 700cabf..f03e88b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 [compat]
 BranchAndPrune = "0.2.1"
 ForwardDiff = "0.10"
-IntervalArithmetic = "0.21"
+IntervalArithmetic = "0.22"
 Reexport = "1"
 StaticArrays = "1"
 julia = "1"
diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index 2e381cc..044ca64 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -34,7 +34,7 @@ const D = derivative
 
 const where_bisect = 0.49609375  # 127//256
 
-const Region = Union{Interval, SVector{N, Interval} where N}
+const Region = Union{Interval, SVector{N, <:Interval} where N}
 
 IntervalBox(x::Interval, N::Integer) = SVector{N}(fill(x, N))
 IntervalBox(xx::Vararg{Interval, N}) where N = SVector{N}(xx...)
diff --git a/src/contractors.jl b/src/contractors.jl
index 3844b00..dd65f07 100644
--- a/src/contractors.jl
+++ b/src/contractors.jl
@@ -46,11 +46,11 @@ function (N::Newton)(X::Interval ; α=where_bisect)
     return m - (N.f(m) / N.f′(X))
 end
 
-function (N::Newton)(X::SVector{<:Interval} ; α=where_bisect)
+function (N::Newton)(X::SVector{M, <:Interval} ; α=where_bisect) where M
     m = interval.(mid.(X, α))
     J = N.f′(X)
     y = gauss_elimination_interval(J, N.f(m))  # J \ f(m)
-    return IntervalBox(m .- y)
+    return m .- y
 end
 
 """
@@ -85,7 +85,7 @@ function (K::Krawczyk)(X::Interval ; α=where_bisect)
     return m - Y*K.f(m) + (1 - Y*K.f′(X)) * (X - m)
 end
 
-function (K::Krawczyk)(X::SVector{<:Interval} ; α=where_bisect)
+function (K::Krawczyk)(X::SVector{N, <:Interval} ; α=where_bisect) where N
     jacobian = K.f′
     mm = mid.(X, α)
     J = jacobian(X)
@@ -100,7 +100,7 @@ end
 Similar to `isempty` function for `IntervalBox`, but also works for `SVector`
 of `Interval`.
 """
-safe_isempty(X) = isempty(IntervalBox(X))
+safe_isempty(X) = any(isempty_interval.(X))
 
 """
     contract(contractor, R)
@@ -113,7 +113,7 @@ function contract(B::Bisection, R::Root)
 
     imX = B.f(X)
 
-    if !(in_interval(0, imX)) || safe_isempty(imX)
+    if !(all(in_interval.(0, imX))) || safe_isempty(imX)
         return Root(X, :empty)
     end
 
@@ -132,7 +132,8 @@ function contract(C::Union{Newton, Krawczyk}, R::Root)
     # Only happens if X is partially out of the domain of f
     safe_isempty(contracted_X) && return Root(X, :unknown)  # force bisection
 
-    NX = intersect_interval(contracted_X, X)
+    NX = intersect_interval(bareinterval(contracted_X), bareinterval(X))
+    NX = IntervalArithmetic._unsafe_interval(NX, min(decoration(contracted_X), decoration(X)), isguaranteed(contracted_X))
 
     !isbounded(X) && return Root(NX, :unknown)  # force bisection
     safe_isempty(NX) && return Root(X, :empty)
@@ -153,7 +154,8 @@ This function assumes that it is already known that `X` contains a unique root.
 """
 function refine(C::Union{Newton, Krawczyk}, X::Region, root_problem)
     while diam(X) > root_problem.abstol
-        NX = intersect_interval(C(X), X)
+        NX = intersect_interval(bareinterval(C(X)), bareinterval(X))
+        NX = IntervalArithmetic._unsafe_interval(NX, min(decoration(C(X)), decoration(X)), isguaranteed(C(X)))
         isequal_interval(NX, X) && break  # reached limit of precision
         X = NX
     end
diff --git a/src/linear_eq.jl b/src/linear_eq.jl
index 22e0f94..f011679 100644
--- a/src/linear_eq.jl
+++ b/src/linear_eq.jl
@@ -143,7 +143,7 @@ function gauss_elimination_interval!(x::AbstractArray, A::AbstractMatrix, b::Abs
         p[i] = (b[i] - temp) / A[i, i]
     end
 
-    return p .∩ x
+    return intersect_interval.(p, x)
 end
 
 function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precondition=true)