diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl
index e3cbb56..b9ae573 100644
--- a/src/IntervalMatrices.jl
+++ b/src/IntervalMatrices.jl
@@ -18,6 +18,7 @@ include("matrix.jl")
 # Arithmetic for interval matrices
 # =================================
 include("arithmetic.jl")
+include("power.jl")
 
 # =======================================================
 # Methods to handle the exponential of an interval matrix
diff --git a/src/power.jl b/src/power.jl
new file mode 100644
index 0000000..089ebf0
--- /dev/null
+++ b/src/power.jl
@@ -0,0 +1,55 @@
+import Base: show
+
+# helper types representing factors
+
+abstract type AbstractFactor{T<:Integer} end
+
+struct NormalFactor{T} <: AbstractFactor{T}
+    a::T
+end
+
+show(io::IO, nf::NormalFactor) = print(io, nf.a)
+
+struct SquaredFactor{T} <: AbstractFactor{T}
+    a::T
+end
+
+show(io::IO, sf::SquaredFactor) = print(io, sf.a, "²")
+
+# square decomposition
+
+function square_decomposition(k::T) where {T<:Integer}
+    factor2source = Dict{T, Vector{AbstractFactor{T}}}()
+    square_decomposition!(factor2source, k)
+    return factor2source
+end
+
+function square_decomposition!(factor2source, k::T) where {T<:Integer}
+    factors = Vector{AbstractFactor{T}}()
+    @assert !haskey(factor2source, k)
+    factor2source[k] = factors
+    a = k
+    while a > 3
+        b = floor(Int, sqrt(a))
+        push!(factors, SquaredFactor(b))
+        if !haskey(factor2source, b)
+            square_decomposition!(factor2source, b)
+        end
+        a -= b^2
+    end
+    if a == 3
+        push!(factors, SquaredFactor(1))
+        if !haskey(factor2source, 2)
+            square_decomposition!(factor2source, 2)
+        end
+        push!(factors, NormalFactor(1))
+    elseif a == 2
+        push!(factors, SquaredFactor(1))
+        if !haskey(factor2source, 2)
+            square_decomposition!(factor2source, 2)
+        end
+    elseif a == 1
+        push!(factors, NormalFactor(1))
+    end
+    return factor2source
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 52ee158..9d41bd6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -80,3 +80,14 @@ end
     f = IntervalMatrices.correction_hull(m, 1e-3, 5)
     f2 = IntervalMatrices.input_correction(m, 1e-3, 5)
 end
+
+@testset "Square decomposition" begin
+    sd = IntervalMatrices.square_decomposition
+    SF = IntervalMatrices.SquaredFactor
+    NF = IntervalMatrices.NormalFactor
+    @test Dict(1 => [NF(1)]) ⊆ sd(1)
+    @test Dict(2 => [SF(1)]) ⊆ sd(2)
+    @test Dict(3 => [SF(1), NF(1)], 2 => [SF(1)]) ⊆ sd(3)
+    @test Dict(9 => [SF(3)], 3 => [SF(1), NF(1)]) ⊆ sd(9)
+    @test Dict(65 => [SF(8), NF(1)], 8 => [SF(2), SF(2)], 2 => [SF(1)]) ⊆ sd(65)
+end