Make behavior of ^(::Matrix{<:Integer},::Integer) consistent with behavior of ^(::Integer, ::Integer)

NEWS.md

@@ -146,6 +146,14 @@ This section lists changes that do not have deprecation warnings.
   * Worker-worker connections are setup lazily for an `:all_to_all` topology. Use keyword
     arg `lazy=false` to force all connections to be setup during a `addprocs` call. ([#22814])
 
+  * `^(A::AbstractMatrix{<:Integer}, p::Integer)` now throws a `DomainError`
+    if `p < 0`, unless `A == one(A)` or `A == -one(A)` (same as for
+    `^(A::Integer, p::Integer)`) ([#23366]).
+
+  * `^(A::AbstractMatrix{<:Integer}, p::Integer)` now promotes the element type in the same
+    way as `^(A::Integer, p::Integer)`. This means, for instance, that `[1 1; 0 1]^big(1)`
+    will return a `Matrix{BigInt}` instead of a `Matrix{Int}` ([#23366]).
+
 Library improvements
 --------------------
 

base/intfuncs.jl

@@ -159,12 +159,19 @@ end
 invmod(n::Integer, m::Integer) = invmod(promote(n,m)...)
 
 # ^ for any x supporting *
-to_power_type(x::Number) = oftype(x*x, x)
-to_power_type(x) = x
-@noinline throw_domerr_powbysq(p) = throw(DomainError(p,
+to_power_type(x) = convert(promote_op(*, typeof(x), typeof(x)), x)
+@noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p,
     string("Cannot raise an integer x to a negative power ", p, '.',
-           "\nMake x a float by adding a zero decimal (e.g., 2.0^$p instead ",
-           "of 2^$p), or write 1/x^$(-p), float(x)^$p, or (x//1)^$p")))
+           "\nConvert input to float.")))
+@noinline throw_domerr_powbysq(::Integer, p) = throw(DomainError(p,
+   string("Cannot raise an integer x to a negative power ", p, '.',
+          "\nMake x a float by adding a zero decimal (e.g., 2.0^$p instead ",
+          "of 2^$p), or write 1/x^$(-p), float(x)^$p, or (x//1)^$p")))
+@noinline throw_domerr_powbysq(::AbstractMatrix, p) = throw(DomainError(p,
+   string("Cannot raise an integer matrix x to a negative power ", p, '.',
+          "\nMake x a float matrix by adding a zero decimal ",
+          "(e.g., [2.0 1.0;1.0 0.0]^$p instead ",
+          "of [2 1;1 0]^$p), or write float(x)^$p or Rational.(x)^$p")))
 function power_by_squaring(x_, p::Integer)
     x = to_power_type(x_)
     if p == 1
@@ -174,9 +181,9 @@ function power_by_squaring(x_, p::Integer)
     elseif p == 2
         return x*x
     elseif p < 0
-        x == 1 && return copy(x)
-        x == -1 && return iseven(p) ? one(x) : copy(x)
-        throw_domerr_powbysq(p)
+        isone(x) && return copy(x)
+        isone(-x) && return iseven(p) ? one(x) : copy(x)
+        throw_domerr_powbysq(x, p)
     end
     t = trailing_zeros(p) + 1
     p >>= t
@@ -196,7 +203,7 @@ function power_by_squaring(x_, p::Integer)
 end
 power_by_squaring(x::Bool, p::Unsigned) = ((p==0) | x)
 function power_by_squaring(x::Bool, p::Integer)
-    p < 0 && !x && throw_domerr_powbysq(p)
+    p < 0 && !x && throw_domerr_powbysq(x, p)
     return (p==0) | x
 end
 

base/linalg/dense.jl

@@ -357,9 +357,15 @@ kron(a::AbstractMatrix, b::AbstractVector) = kron(a, reshape(b, length(b), 1))
 kron(a::AbstractVector, b::AbstractMatrix) = kron(reshape(a, length(a), 1), b)
 
 # Matrix power
-(^)(A::AbstractMatrix, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A), -p) : Base.power_by_squaring(A, p)
+(^)(A::AbstractMatrix, p::Integer) = p < 0 ? power_by_squaring(inv(A), -p) : power_by_squaring(A, p)
+function (^)(A::AbstractMatrix{T}, p::Integer) where T<:Integer
+    # make sure that e.g. [1 1;1 0]^big(3)
+    # gets promotes in a similar way as 2^big(3)
+    TT = promote_op(^, T, typeof(p))
+    return power_by_squaring(convert(AbstractMatrix{TT}, A), p)
+end
 function integerpow(A::AbstractMatrix{T}, p) where T
-    TT = Base.promote_op(^, T, typeof(p))
+    TT = promote_op(^, T, typeof(p))
     return (TT == T ? A : copy!(similar(A, TT), A))^Integer(p)
 end
 function schurpow(A::AbstractMatrix, p)

test/linalg/dense.jl

@@ -582,6 +582,29 @@ end
     end
 end
 
+@testset "issue #23366 (Int Matrix to Int power)" begin
+    @testset "Tests for $elty" for elty in (Int128, Int16, Int32, Int64, Int8,
+                                            UInt128, UInt16, UInt32, UInt64, UInt8,
+                                            BigInt)
+        @test_throws DomainError elty[1 1;1 0]^-2
+        @test (@inferred elty[1 1;1 0]^2) == elty[2 1;1 1]
+        I_ = elty[1 0;0 1]
+        @test I_^-1 == I_
+        if !(elty<:Unsigned)
+            @test (@inferred (-I_)^-1) == -I_
+            @test (@inferred (-I_)^-2) == I_
+        end
+        # make sure that type promotion for ^(::Matrix{<:Integer}, ::Integer)
+        # is analogous to type promotion for ^(::Integer, ::Integer)
+        # e.g. [1 1;1 0]^big(10000) should return Matrix{BigInt}, the same
+        # way as 2^big(10000) returns BigInt
+        for elty2 = (Int64, BigInt)
+            TT = Base.promote_op(^, elty, elty2)
+            @test (@inferred elty[1 1;1 0]^elty2(1))::Matrix{TT} == [1 1;1 0]
+        end
+    end
+end
+
 @testset "Least squares solutions" begin
     a = [ones(20) 1:20 1:20]
     b = reshape(eye(8, 5), 20, 2)

test/linalg/symmetric.jl

@@ -245,10 +245,18 @@ end
                 @testset "pow" begin
                     # Integer power
                     @test (asym)^2   ≈ (Symmetric(asym)^2)::Symmetric
-                    @test (asym)^-2  ≈ (Symmetric(asym)^-2)::Symmetric
+                    if eltya <: Integer && !isone(asym) && !isone(-asym)
+                        @test_throws DomainError (asym)^-2
+                    else
+                        @test (asym)^-2  ≈ (Symmetric(asym)^-2)::Symmetric
+                    end
                     @test (aposs)^2  ≈ (Symmetric(aposs)^2)::Symmetric
                     @test (aherm)^2  ≈ (Hermitian(aherm)^2)::Hermitian
-                    @test (aherm)^-2 ≈ (Hermitian(aherm)^-2)::Hermitian
+                    if eltya <: Integer && !isone(aherm) && !isone(-aherm)
+                        @test_throws DomainError (aherm)^-2
+                    else
+                        @test (aherm)^-2 ≈ (Hermitian(aherm)^-2)::Hermitian
+                    end
                     @test (apos)^2   ≈ (Hermitian(apos)^2)::Hermitian
                     # integer floating point power
                     @test (asym)^2.0   ≈ (Symmetric(asym)^2.0)::Symmetric