From c36a5199f0b8dff0b90118092c8f512750398c78 Mon Sep 17 00:00:00 2001
From: schillic <schillic@informatik.uni-freiburg.de>
Date: Tue, 11 Feb 2020 13:58:08 +0100
Subject: [PATCH 1/2] add square-root decomposition

---
 src/power.jl     | 17 +++++++++++++++++
 test/runtests.jl |  1 +
 2 files changed, 18 insertions(+)

diff --git a/src/power.jl b/src/power.jl
index e62a113..aa8b8a5 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -93,6 +93,7 @@ Increment a matrix power in-place (i.e., storing the result in `pow`).
     * `"multiply"` -- fast computation using `*` from the previous result
     * `"power"` -- recomputation using `^`
     * `"intersect"` -- combination of `"multiply"` and `"power"`
+    * `"sqrt"` -- decompose `k = a² + b`
 
 ### Output
 
@@ -113,6 +114,8 @@ function increment!(pow::IntervalMatrixPower; algorithm::String="intersect")
         pow.Mᵏ = _eval_power(pow)
     elseif algorithm == "intersect"
         pow.Mᵏ = _eval_intersect(pow)
+    elseif algorithm == "sqrt"
+        pow.Mᵏ = _eval_sqrt(pow)
     else
         throw(ArgumentError("algorithm $algorithm is not available; choose " *
             "from 'multiply', 'power', 'intersect'"))
@@ -164,6 +167,20 @@ function _eval_intersect(pow::IntervalMatrixPower)
     return intersect(_eval_multiply(pow), _eval_power(pow))
 end
 
+function _eval_sqrt(pow::IntervalMatrixPower)
+    # decompose k = a² + b with a, b being integers
+    k = pow.k
+    a = floor(Int, sqrt(k))
+    b = k - a^2
+
+    # recursively compute M^a and M^b
+    Mᵏ = square(get(IntervalMatrixPower(pow.M, a)))
+    if b > 0
+        Mᵏ *= get(IntervalMatrixPower(pow.M, b))
+    end
+    return Mᵏ
+end
+
 """
     base(pow::IntervalMatrixPower)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a96c50f..0b59a7d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -142,4 +142,5 @@ end
     increment!(pow)  # power of two
     increment!(pow, algorithm="power")
     increment!(pow, algorithm="intersect")
+    increment!(pow, algorithm="sqrt")
 end

From 6811d03f5074876bcc11374a89c100d052c3357f Mon Sep 17 00:00:00 2001
From: schillic <schillic@informatik.uni-freiburg.de>
Date: Tue, 11 Feb 2020 14:25:20 +0100
Subject: [PATCH 2/2] add algorithm to recursive call

---
 src/power.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index aa8b8a5..96bcca3 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -67,11 +67,12 @@ function IntervalMatrixPower(M::IntervalMatrix{T}) where {T}
     return IntervalMatrixPower(M, M, 1)
 end
 
-function IntervalMatrixPower(M::IntervalMatrix{T}, k::Int) where {T}
+function IntervalMatrixPower(M::IntervalMatrix{T}, k::Int;
+                             algorithm::String="power") where {T}
     @assert k >= 1 "matrix powers must be positive"
     pow = IntervalMatrixPower(M, M, 1)
     @inbounds for i in 1:(k-1)
-        increment!(pow)
+        increment!(pow; algorithm=algorithm)
     end
     return pow
 end
@@ -167,16 +168,16 @@ function _eval_intersect(pow::IntervalMatrixPower)
     return intersect(_eval_multiply(pow), _eval_power(pow))
 end
 
-function _eval_sqrt(pow::IntervalMatrixPower)
+function _eval_sqrt(pow::IntervalMatrixPower; algorithm::String="power")
     # decompose k = a² + b with a, b being integers
     k = pow.k
     a = floor(Int, sqrt(k))
     b = k - a^2
 
     # recursively compute M^a and M^b
-    Mᵏ = square(get(IntervalMatrixPower(pow.M, a)))
+    Mᵏ = square(get(IntervalMatrixPower(pow.M, a; algorithm=algorithm)))
     if b > 0
-        Mᵏ *= get(IntervalMatrixPower(pow.M, b))
+        Mᵏ *= get(IntervalMatrixPower(pow.M, b; algorithm=algorithm))
     end
     return Mᵏ
 end