Matrix power via square decomposition, v2

[schillic]

⚠️ This algorithm is incorrect.

This is an alternative proposal to #83. Apparently pulling the square outside is much better:

julia> A = rand(IntervalMatrix)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-1.28549, 0.82628]    [-0.683239, 0.0622417] 
 [-1.09384, -0.310978]  [-0.162758, 0.00575504]

julia> B = A*A*A*A*A*A*A*A*A
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-121.91, 78.2384]    [-52.3156, 33.4657]
  [-83.7552, 53.5743]  [-35.9422, 22.9383]

julia> C = A^9
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-121.91, 76.1346]    [-52.3156, 32.3239]
  [-83.7552, 43.5264]  [-35.9422, 20.7706]

julia> D = square(square(square(A))) * A  # this is #83
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-121.91, 76.1346]    [-52.3156, 8.98488]
  [-83.7552, 49.3573]  [-35.9422, 20.2105]

julia> sqrt(9)  # this is the proposed approach
3.0

julia> E = square(A^3)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-2.49588, 22.5576]  [-5.84255, 9.67985]
 [-9.35369, 15.4971]  [-2.49588, 6.65116]

julia> diam(B)
2×2 Array{Float64,2}:
 200.148  85.7812
 137.329  58.8804

julia> diam(C)
2×2 Array{Float64,2}:
 198.044  84.6394
 127.282  56.7128

julia> diam(D)
2×2 Array{Float64,2}:
 198.044  61.3004
 133.112  56.1526

julia> diam(E)
2×2 Array{Float64,2}:
 25.0534  15.5224 
 24.8508   9.14703

Of course for non-square powers this is more complicated but can still pay off:

julia> B = A*A*A*A*A*A*A*A*A*A
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-137.338, 213.938]  [-58.9023, 91.8081]
  [-94.296, 146.982]  [-40.3374, 63.0747]

julia> C = A^10
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-113.132, 213.938]   [-53.32, 91.8081]  
  [-84.4426, 146.982]  [-36.1388, 63.0747]

julia> D = square(square(square(A))) * square(A)  # this is #83
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-48.8994, 213.938]  [-52.7449, 91.8081]
 [-84.4426, 146.982]  [-34.1508, 63.0747]

julia> sqrt(10)
3.1622776601683795

julia> E = square(A^3) * A
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-39.5856, 25.0297]  [-16.9877, 2.6562] 
 [-27.1966, 15.5351]  [-11.6708, 6.79703]

julia> diam(B)
2×2 Array{Float64,2}:
 351.275  150.71 
 241.277  103.412

julia> diam(C)
2×2 Array{Float64,2}:
 327.069  145.128 
 231.424   99.2134

julia> diam(D)
2×2 Array{Float64,2}:
 262.837  144.553 
 231.424   97.2254

julia> diam(E)
2×2 Array{Float64,2}:
 64.6152  19.6438
 42.7316  18.4677

[schillic]

This was silly: A^(a²) = A^(a*a) = (A^a)^a != (A^a)² for a != 2.