From 4a68c120d2752843f041095d865cc0972d610c06 Mon Sep 17 00:00:00 2001
From: schillic <schillic@informatik.uni-freiburg.de>
Date: Mon, 27 Jan 2020 15:53:49 +0100
Subject: [PATCH] matrix power

---
 src/IntervalMatrices.jl |  1 +
 src/power.jl            | 86 +++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl        | 23 +++++++++++
 3 files changed, 110 insertions(+)
 create mode 100644 src/power.jl

diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl
index 36da1eb..95db605 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..0e3ffeb
--- /dev/null
+++ b/src/power.jl
@@ -0,0 +1,86 @@
+import Base: ^, show
+
+# helper types representing factors
+
+abstract type AbstractFactor{T<:Integer} end
+
+struct OneFactor{T} <: AbstractFactor{T} end
+
+show(io::IO, ::OneFactor) = print(io, 1)
+
+int(::OneFactor{T}) where {T} = one(T)
+
+struct DoubledFactor{T} <: AbstractFactor{T}
+    a::T
+end
+
+show(io::IO, df::DoubledFactor) = print(io, df.a, "·2")
+
+function int(df::DoubledFactor)
+    return df.a * 2
+end
+
+# binary decomposition
+
+function binary_decomposition(k::T) where {T<:Integer}
+    factor2source = Dict{T, Vector{AbstractFactor{T}}}()
+    binary_decomposition!(factor2source, k)
+    return factor2source
+end
+
+function binary_decomposition!(factor2source, k::T) where {T<:Integer}
+    factors = Vector{AbstractFactor{T}}()
+    @assert !haskey(factor2source, k)
+    factor2source[k] = factors
+    if k == 1
+        push!(factors, OneFactor{T}())
+        return
+    else
+        b = div(k, 2)
+        push!(factors, DoubledFactor(b))
+        if !haskey(factor2source, b)
+            binary_decomposition!(factor2source, b)
+        end
+        if mod(k, 2) == 1
+            # k = 2 * b + 1
+            push!(factors, OneFactor{T}())
+        end
+    end
+    return factor2source
+end
+
+# matrix power
+
+function matpow!(factor2matrix, factor2source, k::Integer)
+    if haskey(factor2matrix, k)
+        return factor2matrix[k]
+    end
+    A = nothing
+    for a in factor2source[k]
+        nₐ = int(a)
+        if !haskey(factor2matrix, nₐ)
+            @assert a isa DoubledFactor
+            # B = C² where C needs to be computed recursively
+            C = matpow!(factor2matrix, factor2source, a.a)
+            B = square(C)
+            @assert !haskey(factor2matrix, nₐ)
+            factor2matrix[nₐ] = B
+        else
+            B = factor2matrix[nₐ]
+        end
+        if A == nothing
+            A = B
+        else
+            A *= B
+        end
+    end
+    factor2matrix[k] = A
+    return A
+end
+
+function ^(M::IntervalMatrix{T1}, k::T2) where {T1, T2<:Integer}
+    factor2source = binary_decomposition(k)
+    factor2source[1] = [OneFactor{T2}()]
+    factor2matrix = Dict{T2, IntervalMatrix{T1}}(1 => M)
+    return matpow!(factor2matrix, factor2source, k)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0e923a7..3968b8e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -91,3 +91,26 @@ end
     b = square(m)
     @test all(inf(a) .<= inf(b)) && all(sup(b) .>= sup(a))
 end
+
+@testset "Binary decomposition" begin
+    bd = IntervalMatrices.binary_decomposition
+    DF = IntervalMatrices.DoubledFactor
+    OF = IntervalMatrices.OneFactor{Int}()
+    @test Dict(1 => [OF]) == bd(1)
+    @test Dict(2 => [DF(1)], 1 => [OF]) == bd(2)
+    @test Dict(3 => [DF(1), OF], 1 => [OF]) == bd(3)
+    @test Dict(9 => [DF(4), OF], 4 => [DF(2)], 2 => [DF(1)], 1 => [OF]) == bd(9)
+    @test Dict(65 => [DF(32), OF], 32 => [DF(16)], 16 => [DF(8)], 8 => [DF(4)],
+               4 => [DF(2)], 2 => [DF(1)], 1 => [OF]) == bd(65)
+end
+
+@testset "Matrix power" begin
+    # flat intervals
+    B = IntervalMatrix([2.0±0 0.0±0; 0.0±0 3.0±0])
+    C = [2.0 0.0; 0.0 3.0]
+    @test all(mid(B^k) == C^k for k in 1:20)
+
+    # proper intervals
+    D = IntervalMatrix([2.0±0.1 0.0±0.1; 0.0±0.1 3.0±0.1])
+    D^4
+end