square decomposition

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

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

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